PaperlayerView exact revision on arXiv

Source-linked AI summary

A second-order linear decoupled scheme with a relaxed zero-energy variable for Rayleigh-B\'enard convection

Matthew McCurdy

arXiv:2610.31364v1physics.flu-dyncs.NAmath.NA

TL;DR

Rayleigh–Bénard convection is continuously driven by heat entering through the lower wall, so energy-stable schemes for dissipative or unheated settings do not address this regime. The paper develops a linear, second-order, fully decoupled scheme with a relaxed zero-energy scalar; its tests show that weighting the scalar equation by the flow’s kinetic energy improves scalar accuracy and permits larger accurate time-steps.

  • Problem

    Existing SAV and GSAV schemes for Boussinesq systems address settings without heat input through the boundary, whereas heating from below supplies energy and limits stability bounds to finite time intervals.

    Source: §1 · pp. 2–3

  • Method

    The paper combines rotational pressure correction, a zero-energy scalar multiplying nonlinear terms, and scalar relaxation in a second-order, linear, fully decoupled scheme for Navier-slip walls.

    Source: §1 · p. 3

  • Results

    Choosing C comparable to kinetic energy E doubled the largest time-step at which steady rolls were reproduced and kept the computational scalar within 5 · 10^-9 of 1 at Ra = 10^6.

    Source: §7 · p. 35

  • Takeaways & Limitations

    Weighting the scalar equation with C comparable to kinetic energy reduces amplification of scalar deviation as flow energy grows, without changing the stability analysis or adding linear solves.

    Source: §7 · p. 35

    The explicit convection still limits the time-step, with δt max |∇u| ≈ 1.2–1.6 in the reported steady-roll tests.

    Source: §6.5 · p. 25

Abstract

from arXiv · show

We construct a second-order, linear, and fully decoupled scheme for Rayleigh--B\'enard convection between Navier-slip walls. The nonlinear terms are treated explicitly and multiplied by a scalar with exact value $1$ and no energy contributed to the continuous problem. As a result, each BDF2 step, with a rotational pressure correction, solves only constant-coefficient linear problems. For the time discretization, we prove that, for any time-step, the energy decays when the walls have equal temperatures, under a condition on the thermal weight, and the energy is bounded on finite time intervals when heated from below. Relaxing the scalar towards $1$ preserves both estimates. Energy stability does not control the accuracy of the scalar; the velocity and temperature fields are computed with its unrelaxed value, whose deviation from $1$ scales like $(E/C)(\delta t/\tau)^3$ for kinetic energy $E$, time-scale $\tau$ of the flow, and weight $C$ of the scalar equation. In our tests, choosing $C$ comparable to $E$ keeps this deviation small and doubles the largest time-step for which steady rolls are reproduced to a relative difference of $10^{-8}$. With Taylor--Hood elements, we observe second-order convergence in time, recover the critical Rayleigh numbers, and match published Nusselt numbers to within $0.5\%$. A three-dimensional implementation passes corresponding validation tests.

1. Introduction

The paper develops a second-order, linear, fully decoupled scheme for Rayleigh–Bénard convection with Navier-slip walls, where heat input prevents global energy decay. It analyzes what the scheme’s zero-energy scalar controls and tests its accuracy in computation.

1. Introduction

  • Unlike prior Boussinesq SAV or GSAV settings without wall heat input, Rayleigh–Bénard convection has a sustained buoyancy energy source and stability bounds only over finite intervals.The paper identifies its setting as heating from below with Navier-slip walls.Source: §1 · p. 3
  • For any time-step, the scheme’s discrete energy decays without heating under a condition on the thermal weight, while heating from below permits a finite-time bound.The analysis concerns the time discretization with exact spatial operators; it does not prove a fully discrete estimate or an error estimate.Source: §1 · p. 4
  • The method combines BDF2, rotational pressure correction, a zero-energy scalar, and scalar relaxation for second-order, linear, fully decoupled solves with Navier-slip walls.Each step reduces to constant-coefficient linear problems; the scalar’s exact value is 1.Source: §1 · pp. 3–4
  • The scalar’s unrelaxed value deviates from 1 at order (E/C)(δt/τ)^3, and computations confirm its dependence on δt^3, 1/C, and kinetic energy E.The practical criterion is C ≳ E; the range of the computational scalar, rather than the stored scalar, diagnoses failure.Source: §1 · p. 4
  • The computational tests include convergence, convection onset, steady-roll heat transport, scalar accuracy and relaxation, energy decay, and a three-dimensional implementation.The tests use Taylor–Hood finite elements; published Nusselt numbers and critical Rayleigh numbers are among the validation targets.Source: §1 · p. 4

2. Model

The model describes Boussinesq convection in a periodic cell between Navier-slip walls, with fixed wall temperatures. Writing temperature as a conductive profile plus a perturbation exposes the energy balance and motivates a zero-energy scalar that preserves the continuous equations while enabling explicit nonlinear terms.

2. Model

  • The Boussinesq model couples incompressible fluid momentum to temperature through buoyancy, with viscosity and thermal diffusion.The Rayleigh and Prandtl numbers specify the numerical experiments’ buoyancy strength and relative momentum and heat diffusion.Source: §2 · p. 6
  • Navier-slip walls impose no penetration and fixed temperatures; the slip coefficient ranges from free slip at zero to no slip in the infinite limit.The domain is periodic in x, with walls at y = 0 and y = H.Source: §2 · p. 6

Conductive steady-state and perturbation.

Energy balance.

  • For heating from below, the buoyancy contribution feeds energy and no positive thermal weight cancels it; with equal wall and reference temperatures, the remaining buoyancy work can be bounded by dissipation.The sign-indefinite energy terms reduce to (G − A T′c)(θ, uy).Source: Energy balance · p. 7

Zero-energy variable.

  • The zero-energy scalar Q multiplies the nonlinear terms, while its continuous equation has zero right-hand side and exact value 1.Multiplying its equation by Q cancels the Q-weighted nonlinear terms in the energy balance; C Q²/2 carries no physical energy, though BDF2 gives it auxiliary numerical energy.Source: Zero-energy variable · p. 7

3. Scheme

The scheme combines BDF2 extrapolation and explicitly scaled nonlinear terms with a rotational pressure correction and a relaxed scalar. The stored scalar supports the energy estimate, while the unrelaxed computational value is used for the velocity and temperature fields.

3. Scheme

  • BDF2 extrapolates velocity and temperature, then computes intermediate velocity, temperature, and scalar values with explicitly treated nonlinear terms.The momentum equation scales the extrapolated convection term by the computational scalar.Source: §3 · p. 8
  • The rotational pressure correction pairs a velocity Helmholtz solve with a Poisson projection and avoids the standard incremental form’s artificial pressure Neumann condition.The Poisson equation uses a homogeneous Neumann condition for the correction potential.Source: §3 · p. 8
  • Only the stored scalar enters the next step; velocity and temperature use the computational scalar, while relaxation moves the stored value toward 1.The relaxation constraint pays for any discrete-energy increase from a fraction η of the step’s dissipation.Source: §3 · pp. 8–9
  • Storing the computational scalar lets explicit-term energy errors accumulate and rescale convection and advection; steady-roll Nu is off by 7% at δt = 5 · 10^-4 with C = 1.The passage attributes the error to the steady state belonging to a different problem.Source: §3 · p. 9
  • The computational scalar’s deviation from 1 scales like (E/C)(δt/τ)^3, and tests confirm the cubic time-step and inverse-weight scaling.This scaling argument is not a bound; C ≳ E is a practical tested rule, not a stability condition.Source: §3 · pp. 9–10
  • At δt ≥ 10^-2, relaxation no longer reset stored Q to 1 with C = 10^4, while it did with C = 100.The passage identifies this regime as beyond the explicit-convection limit.Source: §3 · p. 10

4. Energy stability

This section establishes the scheme’s energy-stability bounds—time-step-independent decay without heating under a thermal-weight condition and a finite-time bound with heating—while clarifying that these results do not establish scalar accuracy or error convergence. It also identifies which parts of the energy argument extend to three dimensions and notes that the flat-wall assumption is needed for Lemma 4.1’s identities.

4.1. Two identities for flat Navier walls

  • Flat Navier-wall identities convert the rotational viscous term into vorticity and wall-friction dissipation, supporting the energy argument for intermediate and divergence-free velocities.The identities rely on flat walls and the Navier-slip condition.Source: §4 · p. 10, §4.1 · p. 11

4.2. Discrete energy and main result

  • The discrete energy includes kinetic, thermal, pressure-correction, and scalar contributions, while the scalar equation cancels the explicit convection and advection terms.The rotational pressure correction contributes a pressure term to the energy.Source: §4 · p. 10, §4.2 · pp. 11, 13
  • For arbitrary time steps, energy decays without heating if Aκ ≥ 10B′, while heating yields a finite-time Grönwall bound for any A > 0; both hold without relaxation.The theorem concerns exact spatial operators.Source: §4 · p. 10, §4.2 · pp. 11–13

4.3. What the theorem does not control

  • The energy bound does not ensure Q stays near 1: relaxation moves stored Q toward 1 but cannot change the unrelaxed value used to compute velocity and temperature.The relaxation contracts the stored value within a step, but this does not guarantee contraction across steps.Source: §4.3 · pp. 14–15

4.4. Why we do not prove an error estimate

  • No error estimate is proved because Grönwall growth, explicit-convection bounds, Navier-slip pressure-correction analysis, and relaxed-scalar error control obstruct it.Numerical evidence instead shows second-order convergence for manufactured solutions and roll formation, and time-step independence of steady states.Source: §4.4 · pp. 15–16

4.5. Three dimensions

  • The two-dimensional identities and energy-stability theorem extend unchanged to a three-dimensional channel.Source: §4.5 · p. 16
  • The dimension-independent analysis—Poincaré inequality in y, BDF2 identities, and pressure terms—extends to three dimensions, with the same constant-coefficient solves per step and one more velocity component.Periodic-face contributions cancel, and the wall-normal flux vanishes because the normal velocity and its tangential derivatives are zero at the wall.Source: §4.5 · p. 16
  • Curved walls add curvature terms to both identities in Lemma 4.1, unlike the flat-wall setting.In two dimensions, the Navier-slip condition becomes ω = (2κ − l/ν) v · τ, and the curvature terms would have to be absorbed by the dissipation.Source: §4.5 · p. 16

5. Implementation

The implementation splits each time step into constant-coefficient field solves, a scalar update, pressure correction, and relaxation. Its finite-element details address energy cancellation and pressure robustness, while the analysis remains limited to exact spatial operators.

5. Implementation

  • The fields use the unrelaxed scalar; after pressure correction, relaxation updates the stored scalar used in the next step.Each step solves four constant-coefficient Helmholtz problems and a pressure Poisson problem before relaxation.Source: §5 · p. 17

Spatial discretization.

  • Taylor–Hood elements use skew-symmetric nonlinear forms, preserving explicit-term cancellation, but the projected velocity is only approximately divergence-free.The velocity and temperature are continuous piecewise quadratic, and pressure is piecewise linear.Source: Spatial discretization · pp. 17–18

Hydrostatic pressure.

  • Removing the conduction-state hydrostatic pressure prevents an O(Gh^2) spurious flow that can exceed the near-onset perturbation being measured.Taylor–Hood elements are not pressure robust, so the unbalanced buoyancy from approximating the quadratic pressure can drive that flow.Source: Hydrostatic pressure · p. 18

Scope of the analysis.

  • The energy analysis covers time discretization with exact spatial operators, not the fully discrete scheme; finite-element energy behavior is checked numerically.The discrete projected velocity does not satisfy the divergence-free condition exactly.Source: Scope of the analysis · p. 18

Reproducibility.

  • A state is deemed steady when bottom-wall Nusselt variation is below 2 · 10^-8 over the run’s final quarter and both wall values agree.Exceptions include free-slip rolls at r = 1.02 and runs with C = 1 and δt = 3 · 10^-3, where residual variation is reported.Source: Reproducibility · p. 18

6. Numerical experiments

The experiments assess convergence, convection onset, steady-roll benchmarks, and the accuracy of the zero-energy variable. With Qn = 1, steady states are independent of the time step and C.

6.1. Manufactured solution

  • The manufactured-solution test shows second-order temporal convergence for velocity, pressure, and temperature; on a wall-resolved mesh, the standard incremental form is measurably less accurate for pressure.On the uniform mesh, the standard form is also second-order, with pressure error about 1.2 times the rotational form.Source: §6.1 · p. 19

Spatial convergence.

  • Spatial orders on the finest meshes approach 3.04 and 3.00 for velocity and temperature in L2, and 2.02, 2.00, and 2.01 for their gradients and pressure.The divergence error decreases at order 2.Source: Spatial convergence · p. 19

6.2. Onset of convection

  • The onset fits give Rac = 657.510 for free-slip walls and Rac = 1707.76 for rigid walls, matching the cited no-slip value 1707.76.The free-slip theoretical value is approximately 657.511.Source: §6.2 · pp. 19, 21

6.3. Steady rolls

  • Near onset, the free-slip steady-roll ratio (Nu − 1)/ϵ decreases from 1.967 to 1.848 as ϵ rises from 0.02 to 0.1, approaching 2.00 as ϵ → 0.The rigid-wall rolls at Ra = 5000 change by only 3 · 10^-6 in Nu between the 32 × 32 and 48 × 48 meshes.Source: §6.3 · p. 22

6.4. Comparison against classic results

  • Present rigid-wall Nusselt numbers are 0.09–0.42% below Clever and Busse's values, while free-slip values agree with Veronis to published rounding through r = 10.The rigid-wall differences are consistent with Clever and Busse's stated error estimates.Source: §6.4 · pp. 23–24, Figure 5 · p. 24

6.5. Accuracy of the zero-energy variable

  • Without relaxation, Q can drift from 1, biasing steady states; at larger steps its collapse can switch off saturating advection and allow unbounded growth.At C = 1, the Nusselt number is off by 7% at δt = 5 · 10^-4.Source: §6.5 · pp. 24–25
  • With C = 100 and relaxation, steady rolls are reproduced to relative difference 10^-8 up to δt = 6 · 10^-3, twice the limit for C = 1.The deviation of ˜Qn from 1 is 100 times smaller and grows like δt^3.Source: §6.5 · p. 25
  • Even with C = 100 or 10^4, explicit convection limits the time step: the observed threshold is δt max |∇u| ≈ 1.2–1.6.Relaxation and C do not remove this restriction.Source: §6.5 · p. 25

Scaling with E

The section examines how the computational scalar and numerical results behave across tested time steps and meshes, including the scalar’s scaling and proximity to 1, code agreement relative to spatial-discretization error, and practical time-step limits despite energy bounds. It also reports validation and three-dimensional tests, while noting that the observed roll sequence is not claimed to be generic.

Scaling with E

  • The fitted exponents were 3.04 for δt, −1.00 for C, and 2.74 for E; the E exponent exceeds the predicted 2.5, with deviations below the scaling line near onset.The fit reproduces runs to a typical factor of 1.25. The E^(5/2) scaling is consistent where H/U is the right flow time scale; C ≳ E is a practical rule, not an analytical condition.Source: Scaling with E · p. 27, Figure 7 · p. 28

6.6. Energy stability versus practical time-step restrictions

  • At Ra = 5000, C = 1 produced a 7% steady-state error at δt = 5·10^-4 and unbounded growth from 10^-3; explicit convection caused oscillations and blow-up at larger steps.With C = 100, ˜Q stayed below 1.2% deviation through δt = 6·10^-3, but the flow-driven explicit-convection limit remained: oscillations appeared near δt = 8·10^-3 and runs blew up at 1.2·10^-2.Source: §6.6 · p. 27

6.7. Energy decay without heating

  • Without differential heating, X^n decreased at every step for δt = 10^-2 to 10^-4, although the largest tested step made peak kinetic energy 17% too large.The test used both A = 5×10^5 and A = 1; the largest step was stable but not accurate.Source: §6.7 · p. 28, Table 10 · p. 28

6.8. Sensitivity to the relaxation parameter

  • For η ≥ 0.5, relaxation reset stored Q^n to 1 each step; with η = 10^-4, Q^n fell to 0.15 at δt = 10^-3 and the solution later blew up.At η = 10^-2 and δt = 3·10^-3, Q^n dipped to 0.939 during the transient, returned to 1, and the steady state was exact. The unrelaxed ˜Q^n is the same for all η ≥ 0.5.Source: §6.8 · p. 29

6.9. A stress test at Ra = 10^6

  • At Ra = 10^6, the successful plume run kept ˜Q^n within 5·10^-9 of 1 and temperature within [0, 1], while the failed run’s ˜Q^n collapsed to 0.32 before blow-up.The successful run used C = 10^6, above the kinetic energy, and δt = 2·10^-6; the authors make no claim about converged statistics, mesh convergence, or behavior beyond t = 0.03.Source: Figure 8 · p. 30, §6.9 · pp. 30–31

6.10. Three dimensions

  • Three-dimensional checks showed second-order temporal convergence, equal growth rates for rolls varying along x or z, and agreement with two-dimensional rolls to the reported precision.The isotropy test gave σ = 3.446650 for both orientations versus 3.446617 in 2D; the 3D and 2D roll Nusselt numbers were 2.1103080 and 2.1103008, respectively.Source: Implementation · p. 31, Table 12 · p. 32

Three-dimensional perturbations of steady rolls

  • A three-dimensional perturbation of steady rolls decayed to at most 3·10^-13 of roll kinetic energy, and the final state agreed with the 2D code to better than 10^-5 relative difference.After t = 0.9, |˜Q^n−1| < 10^-12; the authors describe this as stability to the imposed perturbation in the tested box, not a general stability analysis.Source: Three-dimensional perturbations of… · p. 32, Figure 9 · p. 33

Convection from random initial data

  • In the wide-box run, four straight rolls initially carried 87–93% of the mid-plane spectrum, then gave way to a cellular pattern before diagonal rolls became dominant and steady.Nusselt number settled near 2.10 initially, dropped to 1.85 during the transition, and the final 3D and 2D kinetic energies per unit area were 60.101 and 60.096.Source: Convection from random initial data · p. 33

Convection from random initial data.

  • On a finer 48 × 32 mesh, Nu = 1.950454; comparable meshes differ by 3 · 10−6 in Nu, below the spatial-discretization error of about 5 · 10−5.The stored Q_n is 1 at every step.Source: Convection from random initial data · p. 33
  • Throughout the run, ˜Q_n stays within 1.3 · 10−5 of 1, with its largest deviation at saturation, and within 10−10 once steady.The stored Q_n is 1 at every step.Source: Convection from random initial data · p. 33
  • Random initial data produce straight rolls, a transition, and steady diagonal rolls, but the authors do not claim this pattern sequence is generic.The modes shown are (2, 0) for straight rolls and (1, −1) for the steady diagonal rolls; the box admits only whole-number wavelengths, and stability to other perturbations was not examined.Source: Convection from random initial data · p. 33, Figure 10 · p. 34

7. Concluding remarks

The scheme has finite-time energy bounds for heated convection, but energy stability alone does not ensure accurate scalar values or fields. Tests show that weighting the scalar by kinetic energy improves accuracy, and three-dimensional checks pass in the tested settings.

7. Concluding remarks

  • For any time-step, the time-discrete analysis proves energy decay without heating under a thermal-weight condition and finite-time energy bounds under Rayleigh–Bénard heating.The energy-decay result is not proved for the Taylor–Hood implementation, though its tested unheated runs decreased energy at every step.Source: §7 · p. 35
  • The scheme uses one zero-energy scalar restored by relaxation for Rayleigh–Bénard convection with Navier-slip walls, unlike cited schemes for settings without heat input or with homogeneous boundary conditions.Its per-step solves comprise four constant-coefficient Helmholtz problems and one pressure Poisson problem.Source: §7 · p. 35
  • Energy stability does not ensure scalar accuracy: relaxation corrects stored-variable drift, but the fields use an unrelaxed scalar whose deviation scales with E/C and (δt/τ)^3.Explicit convection also imposes a timestep limit δt max |∇u| = O(1), set by the flow rather than the grid.Source: §7 · p. 35
  • At Ra = 5000, choosing C comparable to E reduced roll-formation error by a factor of two to four and doubled the largest accurate steady-roll timestep.At Ra = 10^6, the same choice kept ˜Q within 5 · 10^-9 of 1 through plume eruption.Source: §7 · p. 35
  • The three-dimensional implementation passes corresponding verification tests; tested roll perturbations decay to round-off, and random-data convection reaches steady diagonal rolls agreeing with the two-dimensional code.In both computations, ˜Q stays within 10^-4 of 1.Source: §7 · p. 35
  • The authors identify larger-box roll instabilities, non-flat or side-heated walls, obstructions, and fluid–porous coupling as possible next steps.The proposed fluid–porous extension is explicitly speculative.Source: §7 · p. 35

Appendix A. The manufactured solution

The manufactured-solution test prescribes a time-dependent exact flow and temperature satisfying the stated wall conditions. Forcing is added to the momentum and heat equations, with the conductive temperature state handled explicitly.

Appendix A. The manufactured solution

  • The test problem uses a periodic channel on [0, 2π) × (−1, 1), with walls at y = ±1 and prescribed material, gravity, and temperature parameters.The walls are held at Tb = 1 and Tt = 0.Source: Appendix A · p. 36
  • The prescribed streamfunction and pressure define an exact solution that is divergence-free and satisfies the no-penetration and wall-temperature conditions.The forcing is added to the momentum and heat equations to realize this solution.Source: Appendix A · p. 36
  • The manufactured Navier-slip condition uses wall forcing gw, with opposite signed expressions prescribed at the lower and upper walls.The supplied expressions are gw = 2(l + ν) sin t sin x at y = −1 and its negative at y = 1.Source: Appendix A · p. 36
  • For the temperature perturbation θ = T − Tc, the heat forcing is fT; the conductive state is time-independent and linear in y, and its buoyancy remains in the momentum equation.The hydrostatic splitting is not used in this test.Source: Appendix A · p. 36

Initial data and steps.

The run starts from zero velocity and pressure with the specified initial temperature, and uses backward Euler for its first step.

Initial data and steps.

Time convergence (Table 2).

The time-error comparison uses a fixed 40 × 40 mesh and L2 differences at t = 0.5 against a finer-step reference run on the same mesh.

Time convergence (Table 2).

Space convergence (Table 3).

The space-error comparison tests four uniform mesh resolutions against exact functions at a fixed time-step.

Space convergence (Table 3).

Post-processing of Figure 8.

Figure 8 is post-processed from computed velocity and temperature fields, with texture indicating flow direction and a bottom panel showing solver-recorded Nusselt numbers.

Post-processing of Figure 8.

Declaration of generative AI and AI-assisted technologies in the manuscript preparation process

The author reports using Claude for language and citation-formatting support, supplemental-code work, and specified technical tasks, while retaining responsibility for the paper.

Declaration of generative AI and AI-assisted technologies in the manuscript preparation process

Funding

The research received no specific grant funding.

Funding

Loading 2610.31364v1…