Bounds on heat transport for convection driven by internal heating
aa r X i v : . [ phy s i c s . f l u - dyn ] F e b This draft was prepared using the LaTeX style file belonging to the Journal of Fluid Mechanics Bounds on heat transport for convectiondriven by internal heating
Ali Arslan † , Giovanni Fantuzzi , John Craske and Andrew Wynn Department of Aeronautics, Imperial College London Department of Civil and Envirnonmental Engineering, Imperial College London(Received xx; revised xx; accepted xx)
The mean vertical heat transport h wT i in convection between isothermal plates drivenby uniform internal heating is investigated by means of rigorous bounds. These areobtained as a function of the Rayleigh number R by constructing feasible solutions toa convex variational problem, derived using a formulation of the classical backgroundmethod in terms of quadratic auxiliary functions. When the fluid’s temperature relative tothe boundaries is allowed to be positive or negative, numerical solution of the variationalproblem shows that best previous bound h wT i / PhysicsLetters A , vol. 377, pp.83-92) can only be improved up to finite R . Indeed, we demonstrateanalytically that h wT i − / R / and therefore prove that h wT i < / R <
65 536.However, if the minimum principle for temperature is invoked, which asserts that internaltemperature is at least as large as the temperature of the isothermal boundaries, thennumerically optimised bounds are strictly smaller than 1 / R = 3 . × .While the computational results suggest that the best bound on h wT i approaches 1 / R → ∞ , we prove that typical analytical constructionscannot be used to prove this conjecture. Key words:
Internally heated convection, turbulent convection, variational methods
1. Introduction
Internally-heated (IH) convection, in which the motion of a fluid is driven by buoyancyforces caused by internal sources of heat, is found in a wide variety of natural andbuilt environments, and plays an essential role in disciplines such as geophysics andastrophysics. For example, radioactive decay drives convection in the Earth’s mantle,which in turn influences plate tectonics and the planet’s magnetic field (Bercovici 2011). Asimilar mechanism explains geological patterns on the surface of Pluto (Trowbridge et al. et al. † Email address for correspondence: [email protected] (Goluskin 2016). IH convection therefore warrants study in its own right to enhancefundamental understanding of buoyancy-driven turbulence, yet has received relativelylittle attention in comparison with Rayleigh–B´enard convection.A fundamental challenge in the study of IH convection, along with many otherturbulent flows, is to characterise the flow’s statistical properties as a function of itscontrol parameters. Following previous work (Goluskin & Spiegel 2012; Goluskin 2016),we consider this problem in an idealized configuration, where a horizontal layer of fluidbetween isothermal plates of equal temperature is heated uniformly at a constant rate.The only control parameters for this setting are the Prandtl number of the fluid, Pr ,and the Rayleigh number based on the internal heating rate, R . Particular statisticalquantities of interest are the dimensionless mean temperature, h T i , and the dimensionlessmean vertical convective heat flux, h wT i , where the mean is obtained by averaging overvolume and infinite time.The dimensionless mean temperature h T i corresponds to the amount of thermaldissipation in the fluid and can be related qualitatively to the proportion of heat withinthe fluid that is transported by conduction, rather than by convection. As described byGoluskin & Spiegel (2012, p.87), the volume integral of outward conduction above andbelow the plane over which the time- and plane-averaged temperature T is maximisedis equal to 2 T at that height. Assuming a homogenised temperature field for sufficientlyhigh R therefore implies that h T i scales in the same way as the maximum T and that theratio of the total (predominantly convective) heat flux to the conductive heat flux, whichcorresponds to a Nusselt number, is 1 / h T i . In contrast, h wT i quantifies the verticalasymmetry caused by the fluid’s motion and is related to the heat fluxes F and F through the top and bottom boundaries by the exact relations (Goluskin 2016) F = 12 + h wT i , F = 12 − h wT i . (1.1)Laboratory experiments (Kulacki & Goldstein 1972; Jahn & Reineke 1974; Mayinger et al. et al. et al. et al. h T i decreases with R . A recent phenomenological theory (Wang et al. h T i ∼ Pr − / R − / when Pr . R − / and h T i ∼ R − / otherwise. The dependence of these predictions on theRayleigh number agrees with rigorous lower bounds on h T i for both finite and infinite P r (Lu et al. a ) up to logarithmic corrections.In contrast to h T i , the behaviour of the mean vertical convective heat flux h wT i remains fascinatingly unclear. While h wT i appears to increase with R in experiments(Kulacki & Goldstein 1972) and three-dimensional simulations (Goluskin & van der Poel2016), it displays little variation with respect to R and non-monotonic behaviour in two-dimensional simulations (Goluskin & Spiegel 2012; Wang et al. h wT i ≈ − σ R − α , (1.2)with prefactors σ ≈ α between 0 .
05 and 0 .
1. Physical theories corrob-orating this power-law behaviour are lacking and the best rigorous mathematical resultavailable to date is the uniform bound h wT i / R -dependent upper bounds on h wT i , some obtainedcomputationally and some proved analytically, the derivation of which relies ontwo key ingredients. The first is a modern interpretation of the background method(Constantin & Doering 1995; Doering & Constantin 1994, 1996) as a particular caseof a broader framework for bounding infinite-time averages (Chernyshenko et al. et al. et al. h wT i even though, contrary to most pastapplications of the background method to convection problems (Doering & Constantin1996; Constantin & Doering 1996; Doering & Constantin 1998; Constantin & Doering1999; Doering & Constantin 2001; Otero et al. et al. et al. b , 2012; Goluskin 2015; Goluskin & Doering2016), this quantity is not directly related to the thermal dissipation. The second keyingredient is a minimum principle, already invoked by Goluskin & Spiegel (2012) andproved in Appendix A, stating that the temperature of the fluid is either no smaller thanthat of the top and bottom plates, or approaches such a state exponentially quickly.Similar results have proved extremely useful in the context of Rayleigh–B´enard convection(Constantin & Doering 1999; Yan 2004; Otto & Seis 2011; Goluskin & Doering 2016;Choffrut et al. R -dependent bounds on h wT i for IH convection at large R .The rest of this work is structured as follows. Section 2 describes the flow configurationand the corresponding governing equations. Heuristic scaling arguments for the meanvertical heat flux are presented in §
3. In §
4, we derive two variational principles to bound h wT i rigorously from above: one that does not consider the minimum principle for thetemperature, and one that enforces it by means of a Lagrange multiplier. Computationaland analytical bounds obtained with the former are presented in §
5, while boundsobtained numerically with the latter are discussed in §
6, along with obstacles to analyticalconstructions. Section 7 offers concluding remarks.Throughout the paper, overlines indicate averages over the horizontal directions andinfinite time, while angled brackets denote averages over the dimensionless volume Ω =[0 , L x ] × [0 , L y ] × [0 ,
1] and infinite time. Precisely, for any scalar-valued function f ( x , t ), f = lim sup τ →∞ τ L x L y ˆ τ ˆ L x ˆ L y f ( x , t ) d x d y d t, (1.3 a ) h f i = lim sup τ →∞ τ ˆ τ Ω f ( x , t ) d x d t, (1.3 b )where ffl Ω ( · )d x = ( L x L y ) − ´ Ω ( · )d x is the spatial average. Note that the f = f ( z )depends only on the vertical coordinate z . We also write k f k and k f k ∞ for the usual L and L ∞ norms of univariate functions f ( z ) on the interval [0 ,
2. Governing equations
We consider a layer of fluid confined between two no-slip plates that are separatedby a vertical distance d and are held at the same constant temperature, which may betaken to be zero without loss of generality. The fluid has density ρ , kinematic viscosity ν , thermal diffusivity κ , thermal expansion coefficient α , specific heat capacity c p , and isheated uniformly at a volumetric rate Q . This corresponds to the configuration denotedby IH1 in Goluskin (2016). To simplify the discussion we assume that the layer is periodicin the horizontal ( x and y ) directions with periods L x d and L y d , but all results presentedin § § L x and L y .To make the problem nondimensional we use d as the length scale, d /κ as the timescale and d Q/κ as the temperature scale. Under the Boussinessq approximation, theNavier–Stokes equations governing the motion of the fluid in the nondimensional domain Ω = [0 , L x ] × [0 , L y ] × [0 ,
1] are (Goluskin 2016) ∇ · u = 0 , (2.1 a ) ∂ u ∂t + u · ∇ u + ∇ p = Pr ( ∇ u + R T ˆz ) , (2.1 b ) ∂T∂t + u · ∇ T = ∇ T + 1 , (2.1 c )with boundary conditions u | z =0 = u | z =1 = 0 , (2.2 a ) T | z =0 = T | z =1 = 0 . (2.2 b )The dimensionless Prandtl and Rayleigh numbers are defined as Pr = νκ , R = gαQd ρc p νκ . (2.3)The former measures the ratio of momentum and heat diffusivity and is a property ofthe fluid, while the latter measures the destabilising effect of internal heating comparedwith the stabilising effect of diffusion and is the control parameter in this study.For any value of R and Pr , equations (2.1a–c) admit the steady solution, u = , T = z (1 − z ), which represents a purely conductive state. This solution is globallyasymptotically stable for any values of the horizontal periods L x and L y when R <
26 926(Goluskin 2016) and is linearly unstable when R is larger than a critical threshold R L ≈
37 325 (Debler 1959), the exact value of which depends on the horizontal periods.Sustained convection ensues in this regime, but has also been observed at subcriticalRayleigh numbers (Tveitereid 1978). Our goal is to characterize the mean verticalconvective heat flux through the layer, h wT i , as a function of R .
3. Heuristic scaling arguments
Phenomenological predictions for the variation of h wT i with the Rayleigh number canbe derived by coupling the total heat budget through the layer with scaling assumptionsfor characteristic length scales δ T and ε T of the lower and upper thermal boundary layers,respectively. These length scales can be defined such that h T i /δ T = F and h T i /ε T = F .Averaging (2.1 c ) over space and infinite time indicates that δ T and ε T satisfy h T i δ T + h T i ε T = 1 , (3.1)while the second identity in (1.1) yields h wT i = 12 − h T i δ T . (3.2)For the sake of definiteness, assume that the mean temperature and δ T decay as powerlaws in R , that is, h T i = R − α /σ and δ T = R − α /σ with α , α >
0. If h wT i approachesa constant as R is raised, then h T i O ( δ T ) and α α , the inequality being strict if h wT i → /
2. Moreover, (3.1) implies that1 ε T = σ R α − σ R α . (3.3)The scalings behind IH convection with the isothermal boundary conditions (2.2) aretherefore necessarily subtle, because the leading scaling of ε T and the correction impliedby δ T both play a crucial role. In particular, any heuristic argument needs to distinguishbetween the physics associated with the unstably stratified flow near the upper boundaryfrom the (very different) stably stratified flow near the lower boundary. In contrast, for IHconvection with an insulating lower boundary, where h wT i = 1 / − T (0), the conjecturedscaling based on empirical observations and lower bounds on h T i leads immediately to T (0) ∼ R − / , corresponding to α = 1 / a ), if the upper boundary layer maintains astate of marginal stability as conjectured for convection by Malkus (1954) and Priestley(1954), then ε T adjusts itself to make the local Rayleigh number Ra ε T (based on theaverage temperature and depth of the upper boundary layer) constant. Expressing Ra ε T in terms of R to leading order, we conclude that Ra ε T = h T i ε T R ∼ σ R − α , (3.4)should be independent of R , which implies that α = 1 / et al. (2020, see regimes III ∞ and IV u ). Alternatively, if oneuses an argument based on balancing a characteristic free-fall velocity p Pr R h T i withthe velocity scale 1 /ε T implied by diffusion at the wall (Spiegel 1963), then to leadingorder ε T p Pr R h T i ∼ σ P r R − α , (3.5)is independent of R , implying that α = 1 /
3. In either case ( α = 1 / α = 1 / ε T anddoes not provide information about the small but crucial correction due to δ T .Conjectures for the scaling of h wT i require an independent argument relating to δ T .Arguably the simplest, but not necessarily the most faithful, comes from assuming thatin some vicinity of the lower boundary there is a balance between heating and diffusionbecause the flow is stably stratified. In terms of the dimensionless variables used here,heating over δ T is proportional to δ T and diffusion is equal to h T i /δ T , which implies that δ T ∼ h T i . This requires α = α /
2, leading to α = 1 / α = 1 / ε T basedon Malkus (1954) or Spiegel (1963), respectively, and therefore to h wT i = 1 / − σ R − /σ or h wT i = 1 / − σ R − /σ . Assuming that max( T ) scales in the same way as h T i ,meaning that the average temperature is approximately uniform away from boundaries,the possibility that δ T ∼ h T i (so α = α /
2) is in good agreement with data fromexperiments and simulations (Goluskin 2016, table 3.2).An alternative argument might consider a Richardson number Ri at the lower bound-ary layer to quantify the destablising effects of turbulence relative to the stabilisingeffects of the density stratification. In terms of dimensionless quantities, the densitystratification is h T i /δ T and the destabilising shear across the lower boundary scalesaccording to p h T i /δ T . Together, these scales imply that Ri ∼ δ T . This is significantbecause, if the flow tends towards a state of marginal stability, then Ri = 1 / h wT i = 1 / − σ R − /σ or h wT i = 1 / − σ R − /σ , corresponding to Malkus (1954) or Spiegel (1963) respectively.The latter scaling would be consistent with the conjectured bound for insulating lowerboundary conditions (Goluskin 2016).
4. Formulation of rigorous bounds
We now turn our attention to bounding h wT i rigorously from above. In §§ § Bounds via auxiliary functional
A rigorous upper bound on h wT i can be derived using the simple observation thatthe infinite-time average of the time derivative of any bounded function V{ u , T } alongsolutions of the governing equations (2.1 a–c ) vanishes, so h wT i = lim sup τ →∞ τ ˆ τ (cid:20) Ω wT + dd t V{ u , T } d x (cid:21) d t. (4.1)In particular, if V can be chosen such that S{ u , T } := U − Ω wT + dd t V{ u , T } d x > a–c ) for some constant U , then h wT i U . The goal,then, is to construct V that satisfies (4.2) with the smallest possible U .While the evolution equations (2.1 b ) and (2.1 c ) cannot be solved explicitly for allpossible initial conditions, when V is given they can be used to derive an explicitexpression for S as a function of u and T alone. Then, to ensure that (4.2) holds alongsolutions (2.1 a–c ) pointwise in time, it suffices to enforce that S{ u , T } be nonnegativewhen viewed as a functional on the space of time-independent incompressible velocityfields and temperature fields that satisfy the boundary conditions, H := { ( u , T ) : (2.2 a,b ), horizontal periodicity and ∇ · u = 0 } . (4.3)We therefore search for a function V and constant U such that S{ u , T } > H . This key relaxation makes this approach tractableand, remarkably, would bring no loss of generality if the governing equations were wellposed (which is not currently known). In this case, in fact, optimizing V over a sufficientlygeneral class of functions would yield an upper bound U exactly equal to the largest valueof h wT i over solutions of (2.1 a–c ) (Rosa & Temam 2020).Unfortunately, the construction of such an optimal V is currently beyond the reach ofboth analytical and computational methods. Nevertheless, progress can be made if werestrict the search to quadratic V in the form V{ u , T } = Ω a Pr R | u | + b | T | − [ ψ ( z ) + z − T d x , (4.4)where the nonnegative scalars a and b and the function ψ ( z ) are to be optimized suchthat (4.2) holds for the smallest possible U . As shown by Chernyshenko (2017), thischoice of V amounts to using the background method (Constantin & Doering 1995;Doering & Constantin 1994, 1996): the profile b [ ψ ( z ) + z −
1] corresponds to the back-ground temperature field, while the scalars a and b are the so-called balance parameters.Note that the z − ψ ( z ), but we isolate it to simplifythe analysis in what follows. Note also that, due to the periodicity in the horizontaldirections, a “symmetry reduction” argument following ideas in Goluskin & Fantuzzi(2019, Appendix A) proves that there is no loss of generality in taking ψ to depend onlyon the vertical coordinate z . Similarly, one can show that the upper bound on h wT i cannot be improved by adding to V a term φ · u , proportional to the velocity field via a(rescaled) incompressible background velocity field φ .To find an expression for S{ u , T } , we calculate the time derivative of the quadratic V in (4.4) using the governing equations (2.1 b,c ), substitute the resulting expressioninto (4.2), and integrate various terms by parts using incompressibility and the boundaryconditions (2.2 a,b ) to arrive at S{ u , T } = Ω (cid:20) aR |∇ u | + b |∇ T | − ( a − ψ ′ ) wT + ( bz − ψ ′ − ∂T∂z + ψ (cid:21) d x + ψ (1) T ′ (1) − ( ψ (0) − T ′ (0) + U − . (4.5)The best bound on h wT i that can be proved with quadratic V is h wT i inf U,ψ ( z ) ,a,b { U : S{ u , T } > ∀ ( u , T ) ∈ H} . (4.6)The right-hand side of (4.6) is a linear optimisation problem because the optimisationvariables, U , ψ ( z ), a and b , enter the constraint S{ u , T } > U linearly.4.2. Fourier expansion
The analysis and numerical implementation of the minimisation problem in (4.6) can beconsiderably simplified by expanding the horizontally-periodic velocity and temperatureas Fourier series, (cid:20) T ( x, y, z ) u ( x, y, z ) (cid:21) = X k (cid:20) ˆ T k ( z )ˆ u k ( z ) (cid:21) e i ( k x x + k y y ) . (4.7)The sums are over wavevectors k = ( k x , k y ) compatible with the horizontal periods L x and L y . We denote the magnitude of each wavevector by k = p k x + k y . The complex-valued Fourier amplitudes satisfy the complex-conjugate relations ˆ u − k = ˆ u ∗ k and ˆ T − k =ˆ T ∗ k , as well as the no-slip and isothermal boundary conditions. Using incompressibilityand writing ˆ w k for the vertical component of ˆ u k , these can be expressed asˆ w k (0) = ˆ w ′ k (0) = ˆ w k (1) = ˆ w ′ k (1) = 0 , (4.8 a )ˆ T k (0) = ˆ T k (1) = 0 . (4.8 b )After inserting the Fourier expansions (4.7) into (4.5) and applying standard estimatesbased on the incompressibility condition and Young’s inequality to replace the horizontalFourier amplitudes ˆ u k and ˆ v k as a function of ˆ w k , the functional S{ u , T } can be estimatedfrom below as S{ u , T } > S { ˆ T } + X k S k { ˆ w k , ˆ T k } , (4.9)where S { ˆ T } = ˆ h b | ˆ T ′ | + ( bz − ψ ′ ) ˆ T ′ + ψ i d z + ψ (1) ˆ T ′ (1) − ( ψ (0) −
1) ˆ T ′ (0)+ U −
12 (4.10)and S k { ˆ w k , ˆ T k } = ˆ (cid:20) aR (cid:18) k | ˆ w ′′ k | + 2 | ˆ w ′ k | + k | ˆ w k | (cid:19) + b | ˆ T ′ k | + bk | ˆ T k | − ( a − ψ ′ ) ˆ w k ˆ T ∗ k i d z. (4.11)Equality holds in (4.9) if u and T depend on only one of the two horizontal coordinates,but not both, and on z .Velocity and temperature fields with a single nonzero Fourier mode are admissible inthe optimization problem (4.6), so the right-hand side of (4.9) is nonnegative if and onlyif each term is nonnegative. Moreover, since the real and imaginary parts of the Fourieramplitudes ˆ w k and ˆ T k give identical and independent contributions to S k , we may assumethem to be real without loss of generality. Thus, we may replace the minimization problemin (4.6) with inf U,ψ ( z ) ,a,b U subject to S { ˆ T } > ∀ ˆ T : (4.8 b ) , S k { ˆ w k , ˆ T k } > ∀ ˆ w k , ˆ T k : (4.8 a,b ) , ∀ k = 0 . (4.12)Any choice of U , ψ ( z ), a and b satisfying the constraints yields a rigorous upper boundon the mean vertical convective heat flux h wT i . Following established terminology, werefer to the inequalities on S k as the spectral constraints.Just like (4.6), the optimization problem in (4.12) is linear and its optimal solutioncan be approximated using efficient numerical algorithms after discretization. For com-putational convenience, however, we simplify the spectral constraints by dropping allnonnegative terms that depend explicitly on k . Specifically, we replace the spectralconstraints with the stronger, but simpler, single condition˜ S{ ˆ w, ˆ T } := ˆ aR | ˆ w ′ | + b | ˆ T ′ | − ( a − ψ ′ ) ˆ w ˆ T d z > ∀ ˆ T , ˆ w : (4.8 a,b ) (4.13)and solve inf U,ψ ( z ) ,a,b U subject to S { ˆ T } > ∀ ˆ T : (4.8 b ) , ˜ S{ ˆ w, ˆ T } > ∀ ˆ w, ˆ T : (4.8 a,b ) (4.14)instead of (4.12). This simplification leads to suboptimal bounds on h wT i at a fixed R but, as discussed in Appendix B, still captures the qualitative behaviour of the optimalones. On the other hand, considering the simplified spectral constraint (4.13) allowsfor significant computational savings when optimizing bounds numerically, because itremoves the need to consider a large set of wavenumbers and it enables implementationusing simple piecewise-linear basis functions (cf. appendix D). This allows for discretiza-tion of (4.14) and its generalization (4.25) derived in § ψ accurately.4.3. Explicit formulation
To simplify the analysis (but not the numerical implementation) of (4.12) it is conve-nient to eliminate the explicit appearance of U . This can be done upon observing that,given any profile ψ ( z ) and balance parameters a , b , the smallest U for which S { ˆ T } isnonnegative over admissible ˆ T is U ∗ = 12 + sup ˆ T k (0)=0ˆ T k (1)=0 (cid:26) − ˆ h b | ˆ T ′ | + ( bz − ψ ′ ) ˆ T ′ + ψ i d z − ψ (1) ˆ T ′ (1) + ( ψ (0) −
1) ˆ T ′ (0) (cid:27) . (4.15)By modifying any admissible ˆ T in infinitesimally thin layers near z = 0 and z = 1, itis possible to show that this value is finite if and only if ψ (0) = 1 and ψ (1) = 0 . (4.15 a,b )Then, the optimal temperature field in (4.15) is given byˆ T ′ ( z ) = ψ ′ ( z ) + 12 b − z , (4.16)and we obtain U ∗ = 12 + 14 b (cid:13)(cid:13)(cid:13)(cid:13) bz − b − ψ ′ ( z ) − (cid:13)(cid:13)(cid:13)(cid:13) − ˆ ψ d z. (4.17)Thus, we may replace the minimization problem (4.12) with the more explicit versioninf ψ ( z ) ,a,b
12 + 14 b (cid:13)(cid:13)(cid:13)(cid:13) bz − b − ψ ′ ( z ) − (cid:13)(cid:13)(cid:13)(cid:13) − ˆ ψ d z subject to ψ (0) = 1 ,ψ (1) = 0 , ˜ S{ ˆ w, ˆ T } > ∀ ˆ w, ˆ T : (4.8 a,b ) . (4.18)Note that although this formulation is not suitable for numerical implementation becausethe cost function is not convex with respect to b , it is more convenient when attemptingto prove an upper bound on h wT i analytically.4.4. Restriction to nonnegative temperature fields
The upper bounding principle derived above can be improved by imposing a mini-mum principle, which guarantees that temperature fields solving the Boussinesq equa-tions (2.1 a–c ) are nonnegative in the domain Ω at large time. More precisely, Appendix Aproves that the fluid’s temperature is nonnegative on Ω at all times if it is so initially,and the negative part of the temperature decays exponentially quickly otherwise.Since h wT i is determined by the long-time behaviour of the velocity and temperaturefields, the minimum principle enables us to replace the upper bound (4.6) with h wT i inf U,ψ ( z ) ,a,b { U : S{ u , T } > ∀ ( u , T ) ∈ H + } , (4.19)where the space H of admissible velocity and temperature fields has been replaced withits subset H + := { ( u , T ) ∈ H : T ( x ) > Ω } . (4.20)The constraint in the modified optimisation problem (4.19) is clearly weaker than theoriginal one in (4.6), so imposing the minimum principle allows for a better boundon h wT i in principle. This is indeed the case at large Rayleigh numbers, as shall bedemonstrated by numerical results in § § S{ u , T } > T is negative on a nonnegligible subset of the domain, we effectively0use a Lagrange multiplier. Specifically, we search for a positive linear functional L suchthat S{ u , T } > L{ T } for all pairs ( u , T ) ∈ H , which satisfy only the boundary conditionand incompressibility. Indeed, the positivity of L implies S{ u , T } > L{ T } > T isnonnegative, as desired. When T is negative on a subset of the domain, instead, L{ T } need not be positive and the constraint on S{ u , T } is relaxed.Analysis in Appendix C proves that there is no loss of generality in taking L{ T } = − Ω q ( z ) ∂T∂z d x , (4.21)where q : (0 , → R is a nondecreasing function, square-integrable but not necessarilycontinuous, to be optimised alongside the bound U , the profile ψ ( z ), and the balanceparameter a, b . If q were differentiable, one could integrate by parts to obtain L{ T } = Ω q ′ ( z ) T d x (4.22)and identify q ′ as a standard Lagrange multiplier for the condition T ( x ) >
0; workingwith (4.21) simply removes the differentiability requirement from q . Moreover, the valueof L{ T } in (4.21) does not change when q is shifted by a constant by virtue of the verticalboundary conditions on T , so we may normalize q such that ˆ q ( z ) d z = ψ (1) − ψ (0) . (4.23)The upper bound (4.19) can therefore be replaced with h wT i inf U,ψ ( z ) ,q ( z ) a,b (cid:26) U : q ( z ) nondecreasing, (4.23) and S{ u , T } + Ω q ( z ) ∂T∂z d x > ∀ ( u , T ) ∈ H , (cid:27) . (4.24)Observe that setting q ( z ) = − ffl Ω q ( z ) ∂T∂z d x to vanish because T iszero at the top and bottom boundaries, so this choice of q results in the upper bound (4.6)derived without the minimum principle for T .The minimisation problem in (4.24) can be expanded using a Fourier series exactly asexplained in § h wT i is bounded above by the optimal valueof the following linear optimization problem:inf U,a,b,ψ ( z ) ,q ( z ) U subject to S { ˆ T } + ˆ q ( z ) ˆ T ′ ( z )d z > ∀ ˆ T : ˆ T (0) = 0 = ˆ T (1) , ˜ S{ ˆ w, ˆ T } > ∀ ˆ w, ˆ T : (4.8 a,b ) ,q ( z ) nondecreasing . (4.25)Furthermore, the constant U can be eliminated from the problem by following the sameprocedure outlined in § S { ˆ T } is minimised whenˆ T ′ ( z ) = ψ ′ ( z ) − q ( z )2 b − z U for given a , b , ψ ( z ) and q ( z ) is U ∗ = 12 + 14 b (cid:13)(cid:13)(cid:13)(cid:13) bz − b − ψ ′ ( z ) + q ( z ) (cid:13)(cid:13)(cid:13)(cid:13) − ˆ ψ d z . (4.27)As one would expect, this expression reduces to (4.17) when q ( z ) = −
1. It is also possibleto show that, since we are only interested in nonnegative ˆ T , the conditions ψ (0) = 1 and ψ (1) = 0 in § ψ (0) ψ (1)
0. Thus, whenthe minimum principle for the temperature is imposed, the optimisation problem (4.18)relaxes into inf ψ ( z ) ,q ( z ) ,a,b
12 + 14 b (cid:13)(cid:13)(cid:13)(cid:13) bz − b − ψ ′ ( z ) + q ( z ) (cid:13)(cid:13)(cid:13)(cid:13) − ˆ ψ d z subject to ψ (0) , ψ (1) ,q ( z ) nondecreasing, (4.23) , ˜ S{ ˆ w, ˆ T } > ∀ ˆ w, ˆ T : (4.8 a,b ) . (4.28)Finally, observe that setting ψ ( z ) = 0, q ( z ) = 0, a = 0 and letting b → / R . Thus, the uniform bound h wT i / V{ u , T } = Ω (1 − z ) T d x , (4.29)which corresponds to the flow’s potential energy measured with respect to the upperboundary. As shown in §
6, however, more general choices can yield better bounds.
5. Optimal bounds for general temperature fields
The best upper bounds on h wT i implied by problem (4.12) can be approximated nu-merically at any fixed Rayleigh number either by deriving and solving the correspondingnonlinear Euler–Lagrange equations (Plasting & Kerswell 2003; Wen et al. et al. R are presentedin § § Numerically optimal bounds
Figure 1 compares the numerically optimal upper bounds U on the mean vertical heattransfer h wT i to the experimental data by Kulacki & Goldstein (1972) and the DNS databy Goluskin & Spiegel (2012) and Goluskin & van der Poel (2016). The bounds werecalculated by solving the minimization problem (4.12) for a fluid layer with horizontalperiods L x = L y = 2, and with the simplified problem (4.14), which is independent ofwavevectors and, therefore, of L x and L y . As expected, the bounds obtained with (4.12)are zero when the Rayleigh number is smaller than the energy stability limit R E ≈
29 723,which differs slightly from the value 26 927 reported by Goluskin (2016) due to ourchoice of horizontal periods. They then increase monotonically with R , showing the samequalitative behaviour as the bounds computed with the simplified problem (4.14), which2 (a)(b) (c) Figure 1. (a)
Optimal bounds on h wT i obtained by solving the wavenumber-dependentproblem (4.12) with horizontal periods L x = L y = 2 ( ) and the simplified problem (4.14)( ). Also plotted are experiments by Kulacki & Goldstein (1972) ( ), 2D DNSs byGoluskin & Spiegel (2012) ( ) and 3D DNSs by Goluskin & van der Poel (2016) ( ). Circlesmark values of R at which optimal profiles ψ ( z ) are plotted in panel (b) . The dashed horizontalline ( ) is the uniform bound h wT i / (b) Optimal profiles ψ ( z ) at R = 10 ( ), 10 ( ) and 10 ( ). (c) Optimal balance parameters a ( , left axis) and b ( , right axis) as afunction of R . reach the value of 1 / R = 259 032. Both sets of results exceed the the uniform upperbound h wT i / ψ ( z ) for the simplified bounding problem (4.14) atselected Rayleigh numbers and the variation of the optimal balance parameters a and b with R are illustrated in Figures 1 (b) and 1 (c) respectively. As expected from thestructure of the indefinite term ( a − ψ ′ ) ˆ w ˆ T of the functional ˜ S , the derivative ψ ′ in the bulkof the domain approaches the value of a as R is raised and leads to the formation of twoboundary layers due to the conditions ψ (0) = 1 and ψ (1) = 0. Moreover, the asymmetryof the boundary layers reflects qualitatively the asymmetry of the IH convection problemwe are studying, which is characterized by a stable thermal stratification near the bottomboundary ( z = 0) and an unstable one near the top ( z = 1). However, note that while ψ is related to the background temperature field, it is not a physical quantity and neednot behave nor scale like the mean temperature in turbulent convection.Other insightful observations can be made by considering the critical temperaturefields ˆ T , which minimize the functional S { ˆ T } for the optimal choice of ψ , a , b and U .These critical temperatures can be recovered upon integrating (4.16), and are plotted inFigure 2 for a selection of Rayleigh numbers. As one might expect, when R is sufficientlysmall such that U = 0, ˆ T = z (1 − z ) coincides with the conductive temperature profile.With the onset of convection and increasing R , boundary layers form at z = 0 and z = 1and the maximum of | ˆ T | decreases. The profiles are also consistent with the uniformrigorous bound h T i /
12 (Goluskin 2016). However, for sufficiently high Rayleigh3 (a) (b) (c)(d) (e) (f )(g)
Figure 2.
Top : Critical temperature ˆ T ( z ) recovered using (4.16). Colors indicate the Rayleighnumber. Panels (a) and (c) show details of the boundary layers. Middle : Detailed view of ˆ T for R = 256 269, 257 500, and 259 032. Dashed lines ( ) are tangent to ˆ T at z = 0. In (d) , ˆ T isnonnegative and has minimum of zero inside the layer. In (e) , ˆ T is initially positive but has anegative minimum. In (f ) , ˆ T ′ (0) = 0 and there is no positive initial layer. Bottom : Upper bounds U on h wT i . Circles mark the values of R considered in (d–f ) and U = 1 / R = 259 032. numbers they are evidently not related to the horizontal and infinite-time averages ofthe physical temperature field, because they become negative near z = 0. Interestingly, asshown in Figure 2 (d) & (e) , this unphysical behaviour first occurs away from the boundaryat R = 256 269, while the numerical upper bound on h wT i reaches the value of 1 / T ′ (0) = 0. The latter is not surprising because the identity T ′ (0) =1 / − h wT i , derived from (1.1) upon recognizing that F = T ′ (0), implies that an upperbound of 1/2 on h wT i is equivalent to a zero lower bound on T ′ (0), which is obtainedwhen ˆ T ′ (0) vanishes. It is therefore clear that the bounding problems (4.12) and (4.14)fail to improve the uniform bound h wT i / R due to a violation of the minimum principle for the temperature, which was nottaken into account when formulating them.5.2. A new Rayleigh-dependent analytical bound
The numerical results in the previous section demonstrate that optimising ψ , a and b cannot improve the uniform bound h wT i at arbitrarily large R . Nevertheless, itis possible to derive better bounds analytically over a finite range of Rayleigh numbersby considering piecewise-linear profiles ψ ( z ) with two boundary layers, such as the onesketched in Figure 3. The numerically optimal profiles in Figure 1, show qualitativelysymmetry with respect to the vertical midpoint z = of the fluid layer, we impose4 zψ ( z )1 1 ψ ′ ( z ) = a δ − εσ Figure 3.
General piecewise-linear ψ ( z ), parametrized by the boundary layer widths δ and ε ,the bulk slope a and the boundary layer height σ . In our proof, we set ε = δ and σ = − a (cid:0) − δ (cid:1) to obtain a profile that is anti-symmetric with respect to z = 1 / anti-symmetry and take ψ ( z ) = − (cid:0) a +12 δ − a (cid:1) z, z δaz + (1 − a ) , δ z − δ (cid:0) a +12 δ − a (cid:1) (1 − z ) , − δ z , (5.1)where δ is a boundary layer width to be specified later. This considerably simplifies thealgebra, at the cost of a quantitatively (but not qualitatively) worse bound on h wT i .Our goal is to determine values for δ , a and b such that ψ satisfies the constraints in thereduced optimisation problem (4.18), while trying to minimise its objective function.To show that the functional ˜ S is nonnegative, observe that the only sign-indefiniteterm in (4.13) is ˆ ( a − ψ ′ ) ˆ w ˆ T d z = a + 12 δ ˆ [0 ,δ ] ∪ [1 − δ, ˆ w ˆ T d z. (5.2)As detailed in Appendix E, this term can be estimated using the fundamental theoremof calculus, the Cauchy–Schwarz inequality, and the boundary conditions (4.8 a,b ) toconclude that ˜ S is nonnegative if δ √ ab ( a + 1) √ R . (5.3)According to the analysis in § a , b and δ satisfying this inequalityproduces the upper bound h wT i
12 + 14 b (cid:13)(cid:13)(cid:13)(cid:13) bz − b − ψ ′ ( z ) − (cid:13)(cid:13)(cid:13)(cid:13) − ˆ ψ ( z ) d z = b
48 + ( a + 1) bδ − ( a + 1) b . (5.4)This bound is clearly minimised when δ is as large as allowed by (5.3), leading to h wT i b
48 + ( a + 1) √ R b √ ab − ( a + 1) b . (5.5)5 Figure 4.
Comparison between the semi-analytical bounds computed with (5.5) and (5.6 a,b )( ), the explicit bound (5.8) ( ), the numerically optimal bounds obtained in § L x = L y = 2) and with (4.14) ( , • ), and experimental and DNS data (see Figure 1 for akey of the symbols). A dashed vertical line ( ) indicates the smallest energy stability threshold, R ≈
26 926, while a dashed horizontal line ( ) indicates the uniform upper bound h wT i . The values of a and b minimising the right-hand side of this inequality satisfy −√ b + √ R √ a ( a + 1)(5 a −
1) = 0 , (5.6 a )1 + 12( a + 1) b − a + 1) √ R √ ab = 0 , (5.6 b )and can be computed numerically for fixed R to obtain upper bounds plotted as a dot-dashed line in Figure 4.Slightly worse but fully analytical bounds can instead be obtained if we drop the lastterm from (5.5). In this case, the optimal a and b are found explicitly as a = 15 and b = 95 (cid:18) R (cid:19) (5.7)such that we obtain h wT i − R . (5.8)This bound, plotted as a solid line in Figure 4, is smaller than the uniform bound of1/2 up to R = 2 = 65 536, which is approximately 2.43 times larger than the energystability threshold.
6. Optimal bounds for nonnegative temperature fields
The analysis and computations in § h wT i / R , one must invoke the minimum principle forthe temperature explicitly to avoid unphysical critical temperatures ˆ T that are negativein the interior of the layer. As discussed in § h wT i that take thisconstraint into account can be found by solving (4.25).6.1. Numerically optimal bounds
Problem (4.25) was discretised into an SDP and solved with the high-precision solver sdpa-gmp (Yamashita et al. . × R . × . The MATLAB toolbox SparseCoLO (Fujisawa et al. (a) (b) (c)
Figure 5. (a)
Optimal bounds U on h wT i computed by solving (4.25), which incorporatesthe minimum principle for temperature ( ). Also plotted ( ) are the bounds computed in § q ( z ) = − (b) Detail of the region inside the red dashed box( ) in panel (a) . (c) Difference between our optimal bounds and the uniform bound h wT i / (b) and (c) , a circle at R ≈
256 269 marks the pointat which q ( z ) begins to vary from −
1, while stars ( ∗ ) mark the Rayleigh numbers at which ψ are plotted in Figure 6. Rayleigh number, we employed the finite-element discretisation approach described inAppendix D on a Chebyshev mesh with at least 6000 piecewise-linear elements, increasingthe resolution until the upper bounds changed by less than 1%. Achieving this at R =3 . × required approximately 12 200 elements. The numerical challenges associatedwith setting up the SDPs accurately in double-precision using SparseCoLO on evenfiner meshes prevented us from considering a wider the range of Rayleigh numbers.Figure 5 compares numerical upper bounds on the heat transfer obtained with (solidline) and without (dot-dashed line) the minimum principle for the temperature, that is,by optimising q ( z ) or by setting q ( z ) = − § q ( z ) = − R <
256 269. For higher R , the numerical upperbounds on h wT i with optimised q ( z ) are strictly better than those with q ( z ) = − from below as R is raised. Notethat although the deviation from is small at the highest values of R that could behandled, it is much larger than the tolerance (10 − ) used by the multiple-precision SDPsolver sdpa-gmp , giving us confidence that it is not a numerical artefact. Moreover, thedeviation from of the numerical bounds, shown in Figure 5 (c) , appears to decay as apower law, suggesting that the optimal bound available within our bounding frameworkmay have the functional form (1.2). Although the range of Rayleigh numbers spannedby our computations is too small to accurately predict the exponent and the prefactor,it is clear that the decay is much faster than predicted by the heuristic arguments in § ψ ( z ) for selected Rayleigh numbers are shown in Figures 6 (a,b) anddiffer significantly from the corresponding profiles in Figure 1 (b) obtained when theminimum principle for the temperature is disregarded. When the minimum principleis enforced, ψ appears to approach zero almost everywhere as R is raised, but alwayssatisfies ψ (0) = 1. This leads to the formation of a very thin boundary layer near z = 0,which at high R consists of three distinct sublayers identified by two points, z = δ and z = δ , at which ψ is not differentiable. These are indicated by gray dashed lines inFigure 6 (b) . In the first sublayer, from z = 0 to z = δ , ψ is observed to vary linearly.The second sublayer, δ < z < δ , is observed only for R > . × and we observethat ψ ∼ z − approximately. The third sublayer, from z = δ to the point z = δ at7 (a) (b)(c) (d) Figure 6. (a)
Optimal ψ ( z ) in the entire domain plotted from 0 ψ . ψ (0) = 1for all R . (b) Logarithmic plot to highlight behaviour near zero with 0 ψ
1, the 4 colourplots correspond to the R highlighted in Figure 5, for R = 2 . . . . × . The dashed lines ( ) highlight the boundaries of the sub layers (0 , δ ) and ( δ , δ ). (c) Variation with R of the optimal balance parameters a ( ) and b ( ) for (4.25) and (d) theratio of balance parameters b/a when the Lagrange multiplier is active. which ψ attains a local minimum, does not have a simple functional form. In the bulk, ψ increases approximately linearly with slope very close to a , as was the case in §
5, andthe condition ψ (1) = 0, which emerges as a result of the optimization and is not imposed a priori , is attained through a small boundary layer of width ε near z = 1. We choosethe boundary of this layer as the point z = 1 − ε at which ψ has a local maximum.Figure 6 (c) shows the variation of the balance parameters with R . Below 256 269, bothcoincide with the values plotted in Figure 1 (c) . At higher R , the minimum principle forthe temperature becomes active and both balance parameters start to decay rapidly. Itwould be tempting to conjecture that a ∼ b ∼ R − p for some power p but, given the smallvariation of a/b evident in Figure 6 (d) , we cannot currently exclude that subtly differentscaling exponents or higher-order corrections do not play an important role in obtainingan upper bound on h wT i that approaches asymptotically from below.Figure 7 shows the structure of q , whose (distributional) derivative represents theLagrange multiplier enforcing the minimum principle. The multiplier, therefore, is activein regions where q is not constant. The choice q ( z ) = − R = 256 269.At higher Rayleigh numbers, the multiplier becomes active between z = δ and z = δ ,which is what causes the second boundary sublayer in ψ . In the immediate vicinity ofthe bottom boundary (0 z δ ), q is constant and it appears that q ( z ) − ψ ′ ( z ) ≈ b/ (c) ). Indeed, inspection of the cost function in (4.28) suggests that q ( z ) − ψ ′ ( z ) = b (1 − z ) / z -dependent correction in our numerical results. In the second sublayer ( δ z δ ),where the Lagrange multiplier is active, q ( z ) ∼ − z − . Again, this is consistent with the8 R . × . × . × . × . × Figure 7.
Variation of q ( z ) shown on a logarithmic scale to highlight the non-zero region ofthe Lagrange multiplier q ′ ( z ), from R = 2 . − . × . The colorbar here applies for all ψ , q and ˆ T presented in § (a) (b) Figure 8.
Variation with R of the critical temperature profiles, ˆ T , plotted in the region of0 z . × − , (a) plots of individual ˆ T from R = 2 . − . × , (b) contour plotdepicting ˆ T as R increases, the solid white line corresponds to the region within which q ′ ( z ) > T = 0, the colorbar illustrates the value of ˆ T , values below 10 − are to be assumed zero. minimization of the cost function in (4.28), as one expects q to cancel the very largecontribution of ψ ′ near the bottom boundary.Further validation of our numerical results comes from inspection of the criticaltemperatures ˆ T ( z ), which can be recovered using (4.26) and are shown in Figure 8 (a,b) for selected values of R . As expected, the critical temperatures are nonnegative for all z and vanish identically (up to small numerical tolerances) in the region δ z δ , where q is active. We note that, for a given Rayleigh number, this region is strictly larger thanthe range of z values for which the critical temperatures in Figure 2 (c) are negative,indicating that the minimum principle alters the problem in a more subtle way thansimply saturating the constraint T > χ ( z ) = ψ ′ ( z ) − q ( z ) , (6.1)and rewrite (4.27) as h wT i
12 + b " ˆ (cid:18) z − − χ ( z ) b (cid:19) − ψ ( z ) b d z . (6.2)Panels (a) , (c) & (d) in Figure 9 suggest that χ (0) /b → − and χ (1) /b → − as R → ∞ .In fact, profiles χ ( z ) /b for different Rayleigh numbers collapse almost exactly throughout9 (a) (b)(c)(d) (e) Figure 9. (a)
Profiles χ ( z ) /b for 2 . × R . × . (b) Detailed view of χ ( z ) /b inthe bulk. (c,d) Variation of χ (0) /b and χ (1) /b with R . (e) Contributions to the integral of ψ/b in the regions (0 , δ ) ( ), ( δ , δ ) ( ), ( δ , δ ) ( ), ( δ , − ε ) ( ) and (1 − ε,
1) ( ), as afunction of R . the layer, but subtle corrections are present; for instance, Figure 9 (b) demonstrates that,in the bulk, the mean value of χ ( z ) /b decreases with R . It is also evident that χ ( z ) is notexactly constant throughout the bulk, but increases by approximately 10 − . Since q ( z )is constant in this region, we conclude that ψ ′ ( z ) is not constant, but displays subtle andthus nontrivial variation.Figure 9 (e) illustrates the variation with R of contributions to the integral of ψ ( z ) /b from regions (0 , δ ), ( δ , δ ), ( δ , δ ), ( δ , − ε ) and (1 − ε, δ z − ε ), but it slowly decreases with R . The same is true ofthe contribution of the top boundary layer (1 − ε z
1) and the outermost boundarysublayer near z = 0 ( δ z δ ). Only in the first two boundary layers near z = 0 doesthe value of the integral increase with R , suggesting that the integral of ψ ( z ) /b near theboundary layer may become the dominant term as R → ∞ . While the range of Rayleighnumbers covered by our computations is too small to confirm or disprove this conjecture,it is certain that the integral of ψ ( z ) /b must remain large enough to offset the positiveterm in (6.2) in order to obtain a bound on h wT i smaller than 1 / R , q ( z ) = ψ ′ ( z ) + b , z δ , − q z − , δ z δ ,ψ ′ (1) + b , δ z , (6.3)for some positive constant q , and that a, ψ, q ∼ b . The next section investigates whethersuch ansatz can lead to upper bounds on h wT i that are strictly smaller than 1 / Towards an analytical bound
The numerical evidence presented in § h wT i obtained from the optimization problem (4.28) approaches asymptotically from belowas R → ∞ . To confirm this observation with a proof requires, for every R > a, b, ψ ( z ) , q ( z ) whose corresponding cost isstrictly less than . This section discusses the challenges presented by this goal.To aid the discussion, the following definition introduces three subsets A , B and C ofdecision variables which capture some of their properties observed from the numericalstudy in § Definition 1.
Let < δ < − ε < and γ > and R > . Decision variables ( a, b, ψ ( z ) , q ( z )) are said to belong to:1. The set A{ δ, ε } if the following conditions hold:(a) The balance parameters satisfy < a b ;(b) ψ ∈ C [0 , with boundary conditions ψ (0) and ψ (1) = 0 ;(c) The derivative of ψ satisfies ψ ′ | [0 ,δ ] , ψ ′ | [ δ, − ε ] > , and k ψ ′ k L ∞ (1 − ε, b .2. The set B{ γ } if both ψ ′ ( z ) and q ( z ) are constant on the interval ( − γ, + γ ) .
3. The set C{ δ, ε, R } if δ k a − ψ ′ k L ∞ (0 ,δ ) + ε k a − ψ ′ k L ∞ (1 − ε, r abR . The set A{ δ, ε } contains profiles ψ which possess an initial (and potentially severe)boundary layer in an interval [0 , δ ], then increase in a bulk region [ δ, − ε ], beforeapproaching ψ (1) = 0 in an upper boundary layer contained in the interval [1 − ε,
1] andin which ψ ′ is controlled by the balance parameter b , as seen in previous results. Optimaldecision variables ( a, b, ψ ( z ) , q ( z )) obtained in § A{ δ, ε } with the exception of the differentiability condition ψ ∈ C [0 , ψ appear to be piecewise differentiable, losingdifferentiability at two points corresponding to the boundaries of non-constant behaviourof the multiplier q ( z ) observed in figure 7. However, since no higher derivatives of ψ appearin the optimization problem (4.28), adding the constraint ( a, b, ψ ( z ) , q ( z )) ∈ A to (4.28)will not change its optimal cost.The set B{ γ } contains decision variables for which ψ ′ and q are constant in someinterval centred at . Figures 7 and 9 reveal that this is not the case for the numericallyoptimal ψ , so the use of B{ γ } corresponds to a proof which ignores subtle variationsfrom a purely linear profile away from the boundaries. Without further assumptions (say,restriction to fluids with infinite Prandtl number), it is not clear how such variations canbe exploited in analytical constructions.The set C{ δ, ε, R } relates to a choice of profile ψ and balance parameters a, b for whichthe constraint ˜ S > R . In particular, if a, b, ψ ( z ) satisfy(E 8) and it is the case that ψ ′ | [ δ, − ε ] = a , then the argument in Appendix E implies that˜ S >
0. Specifically, C{ δ, ε, R } provides sufficient control of the severity of the boundarylayers of ψ for the spectral constraint to be provably satisfied. While crude, estimates ofthis form in conjunction with constant ψ ′ ( z ) in a bulk region δ z − ε are employedfor almost all analytical constructions of background fields.We now return to the original question of attempting to upper-bound h wT i via ananalytical construction of feasible decision variables ( a, b, ψ ( z ) , q ( z )) for (4.28). It isnot unreasonable, based upon the above evidence, to propose R -dependent balanceparameters a = a R , b = b R , boundary layer widths δ = δ R , ε = ε R and profiles ψ R , q R a R , b R , ψ R ( z ) , q R ( z )) ∈ A{ δ R , ε R } ∩ B{ γ } ∩ C{ δ R , ε R , R } , R > , for some γ >
0. The following result shows that, for such a construction, there is a hardlower bound on the optimal cost achievable using (4.28).
Proposition 1.
Let < δ < − ε < and γ > , R > . Suppose that ( a, b, ψ ( z ) , q ( z )) ∈ A{ δ, ε } ∩ B{ γ } ∩ C{ δ, ε, R } . (6.4) Then b (cid:13)(cid:13)(cid:13)(cid:13) bz − b − ψ ′ ( z ) + q ( z ) (cid:13)(cid:13)(cid:13)(cid:13) − ˆ ψ d z > b γ − · √ R ! . From the proof in Appendix F, the consequence of Proposition 1 is the following. If oneconstructs feasible decision variables for the optimization problem (4.28) which satisfy(6.4), then the best achievable bound U for which h wT i U must satisfy U >
12 + b γ − · √ R ! . Hence, using such a construction with the assumption that ψ ′ and q are constant in abulk region ( − γ, + γ ), it is not possible to prove that h wT i < at arbitrarily highRayleigh number.Consequently, one must ask what conditions should be dropped if a rigorous bound h wT i < , R >
1, is to be found. The numerical evidence presented in § A . Consequently, either B or C must bedropped. Figure 7 indicates that optimal multipliers q ( z ) are constant outside a lowerboundary layer. Hence, dropping B corresponds to choosing a profile with non-constant ψ ′ ( z ) in the bulk; the cost function of (4.28) and Figure 9 (b) suggests that a quadraticansatz for ψ ( z ) may be beneficial. Dropping C corresponds to requiring more sophisticatedanalysis of the spectral constraint. Using these insights will be the focus of future research.
7. Conclusions
We obtain upper bounds U on the vertical heat transport h wT i in IH convection thatimprove on the best existing bound h wT i / R . Specifically, we construct an analytical proof that h wT i − R / ,which improves on the uniform bound for R < = 65 536 (2.43 times larger than theenergy stability threshold). However, our numerical results suggest that the best availablebound tends to 1 / R → ∞ when a Lagrange multiplier isintroduced to enforce non-negativity of the optimal temperature field.A numerical challenge encountered in this work was the implementation of the optimi-sation problem (4.25). The sharp boundary layers in ψ and q near z = 0 require extremelyhigh resolution and, consequently, a fine mesh. This poses challenges not only in termsof computational cost, but also in terms of numerical accuracy. Using a simplificationof the spectral constraint we obtained results for a limited range of R , but in order toshed light on the problems highlighted above, possible improvements to the algorithm orproblem’s formulation remain a point of interest and opportunity for further research.2As explained in § U → / R → ∞ requires the use of more sophisticated estimates to satisfy thespectral constraint than those that are typically used. Such improvements are necessary,rather than sufficient, conditions because there is a possibility that U = 1 / R that lie outside the range that we have studied numerically. In this regard, thenumerical results, which correspond to the best available bound for quadratic auxiliaryfunctions, suggest either exponential or extremely strong algebraic scaling of U towards1 / R that is far from any of the heuristic scaling possibilities discussed in §
3. Moresophisticated treatments of the spectral constraint might yield a proof that U → / § et al. et al. Acknowledgements
We are grateful to D. Goluskin for enlighting us on many aspects ofIH convection and for sharing his DNS data. This work also benefited from conversationswith C. Nobili, J. Whitehead, C. R. Doering, I. Tobasco and G. O. Hughes. GF gratefullyacknowledges the support of an Imperial College Research Fellowship and the hospitalityof the 2018 GFD program at Woods Hole Oceanographic Institution, where this workwas begun.
Conflict of interests
The authors report no conflict of interests.
Appendix A. Minimum principle in IH convection
A minimum principle for IH convection can be proved by adapting arguments forRayleigh–B´enard convection (Foias et al. T − ( x , t ) := max {− T ( x , t ) , } (A 1)denote the negative part of T and observe that T − is nonnegative on the fluid’s domain Ω . Multiplying the advection-diffusion equation (2.1 c ) by T − and integrating by partsover the domain using the boundary conditions and incompressibility yields12 dd t k T − k = −k∇ T − k − ˆ Ω T − d x . (A 2)Upon observing that the last integral on the right-hand side is positive and that T − vanishes at the top and bottom boundaries, so k∇ T − k > µ k T − k for some constant µ > t k T − k − µ k T − k . (A 3)Gronwall’s lemma then yields k T − ( t ) k k T − (0) k e − µt , (A 4)so T − tends to zero in L ( Ω ) at least exponentially quickly. This implies that T ( x , t ) > T ( x , t ) is nonnegative at all times ifit is so at t = 0.3 (a) (b) (c)(d) (e) (f )(g) Figure 10.
Top : Critical temperature ˆ T ( z ) recovered using (4.16) when the spectral constraintis given by S k >
0. Colors indicate the Rayleigh number. Panels (a) and (c) show details ofthe boundary layers.
Middle : Detailed view of ˆ T for R = 7 353 761, 7 455 000, and 7 600 269.Dashed lines ( ) are tangent to ˆ T at z = 0. In (d) , ˆ T is nonnegative and has minimum of zeroinside the layer. In (e) , ˆ T is initially positive but has a negative minimum. In (f ) , ˆ T ′ (0) = 0and there is no positive initial layer. Bottom : Upper bounds U on h wT i . Circles mark the valuesof R considered in (d–f ) and U = 1 / R = 7 600 269. Appendix B. Comparing the full and simplified spectral constraints
This appendix provides further computational evidence that replacing the spectral con-straints S k > S > h wT i . The simplified optimization problem (4.14)was solved using the finite-element expansion approach described in appendix (D). Forsimplicity, instead, the wavenumber-dependent problem (4.12) was solved using thegeneral-purpose toolbox quinopt (Fantuzzi et al. Lemma 1.
Fix R, a , b and ψ ( z ) satisfying ψ (0) = 1 and ψ (1) = 0 . The inequality S k { ˆ w k , ˆ T k } > holds for all ˆ w k and ˆ T k satisfying conditions (4.8a,b) if k > R ab k a − ψ ′ k ∞ . (B 1) Proof.
The H¨older and Cauchy–Schwarz inequalities yield (cid:12)(cid:12)(cid:12)(cid:12) ˆ ( a − ψ ′ ) ˆ w k ˆ T k d z (cid:12)(cid:12)(cid:12)(cid:12) k a − ψ ′ k ∞ k ˆ w k k (cid:13)(cid:13) ˆ T k (cid:13)(cid:13) . (B 2)4 (a) (b) Figure 11. ( a ) Optimal ψ ( z ) for selected Rayleigh numbers (see colorbar) obtained by solvingoptimization problem (4.12), which imposes the exact spectral constraints S k >
0. ( b ) Optimal ψ ( z ) at the same Rayleigh numbers but obtained by solving optimization problem (4.14), whichimposes the simplified spectral constraint ˜ S > Consequently, S k > ak R k ˆ w k k − k a − ψ ′ k ∞ k ˆ w k k (cid:13)(cid:13) ˆ T k (cid:13)(cid:13) + bk (cid:13)(cid:13) ˆ T k (cid:13)(cid:13) . (B 3)The right-hand side is a homogeneous quadratic form in k ˆ w k k and k ˆ T k k and isnonnegative for any choice of ˆ w k and ˆ T k if and only if (B 1) holds.This result guarantees that, when implementing (4.12) numerically, it suffices toconsider wavenumbers with k < (cid:18) R ab (cid:19) k a − ψ ′ k ∞ . (B 4)While the right-hand side of this inequality is unknown a priori, as it depends on theoptimisation variables a , b and ψ , in practice one can simply solve (4.12) using allwavevectors with k smaller than an arbitrarily chosen value. Then, one checks whether S k is indeed nonnegative for all k satisfying (B 4), and repeats the computation witha larger set of wavevectors if these checks fail. Upper bounds obtained by solving thefull problem (4.12) and the simplified problem (4.14) are shown in Figure 1. Here wedemonstrate the equivalence of the results for both spectral constraints qualitatively. Asexpected, using the simplified spectral constraint yields worse bounds at a fixed Rayleighnumber. While the full spectral constraint yields bounds that are zero up for all R upto the energy stability threshold, which depends on the choice of horizontal periods L x and L y , the simplified functional ˜ S is insensitive to these values and gives a conservativeestimate for the nonlinear stability threshold. Nevertheless, both sets of result displaythe same qualitative increase as R is raised. This is further demonstrated in Figures 10,which can be compared to Figure 2In particular, the upper bound reaches 1 / T ,which minimizes the functional S , has zero gradient at z = 0, as can be observed inFigures 10 (f ) & (g) . Shown in panels (a)-(c) are the critical temperature fields, ˆ T , for10 R . In the middle row of the figure, going from left to right, observe firstthat for R > (e) shows a R at which ˆ T is positive in a small region very close to the wall but clearlyviolates the minimum principle further away. In (f ) , where R = 7 600 269, for our choice5of L x = L y = 2, the numerically optimal bound passes the uniform bound, at which pointwe have ˆ T ′ (0) = 0. These results for the k -dependent spectral constraint qualitativelymatch the results for the simplified spectral constraint (4.13) presented Figure 2.As evidenced by Figure 11, the optimal ψ obtained with the full and simplified spectralconstraints also display similar features. The only notable difference is the non simplebehaviour near the boundary layer of the optimal profiles ψ for the full problem (4.12),which arise after bifurcations in the critical wavenumbers as R is raised. In panel ( b ),instead, ψ exhibits a simple structure in both of the boundary layers.These observations confirm that strengthening the spectral constraints using thewavevector-independent functional ˜ S only affects our computational results quantita-tively, but preserves the overall qualitative behaviour. Appendix C. Justification of (4.21)
To justify the choice of linear functional in (4.21), we start with a technical lemma.In this appendix, T is the space of square-integrable temperature fields in the Sobolevspace H ( Ω ) that are horizontally periodic and vanish at z = 0 and 1. We equip T withthe inner product ( T , T ) = ffl Ω ∇ T · ∇ T d x . Lemma 2.
Let q : (0 , → R be a square-integrable function. The linear functional L q ( T ) = − Ω q ( z ) ∂T∂z d x on T is positive if and only if q ( z ) is nondecreasing.Proof. First, we prove that q is nondecreasing if L q is a positive functional. Fix any z , z ∈ (0 ,
1) with z < z and choose ε > < z − ε and z + ε < T ε ( x ) = T ε ( z ) that varies only in z and satisfies ∂T ε ∂z := ε − z − ε z z + ε, − ε − z − ε z z + ε, . Clearly, T ε ∈ T and is nonnegative, so the positivity of L q yields0 L q { T ε } = − Ω q ( z ) ∂T ε ∂z d x = 1 ε ˆ z + εz − ε q ( z )d z − ε ˆ z + εz − ε q ( z )d z. (C 1)Letting ε → q ( z ) q ( z ), which implies that q is nondecreasing since z and z are arbitrary.To prove the reverse statement, suppose that q is nondecreasing but that L q is notpositive. This means that there exist a constant c > T ∈ T ,nonnegative on the domain Ω , such that L q ( T ) − c . By a standard approximationargument, we may also assume that q ( z ) is smooth on [0 , T yields Ω q ′ ( z ) T d x = L q ( T ) − c. This is a contradiction because the left-hand side is a nonnegative quantity, as q ′ ( z ) > q is nondecreasing) and T ( x ) > Ω by assumption.6Lemma 2 guarantees that the linear functional in (4.21) is positive, as required, if q is nondecreasing. The next result shows that considering more general types of positivelinear functionals is not helpful. Proposition 2.
Suppose there exists a positive linear functional L : T → R such that S{ u , T } > L{ T } for all pairs ( u , T ) ∈ H . Then, there exists a nondecreasing square-integrable function q : (0 , → R such that S{ u , T } > − Ω q ( z ) ∂T∂z d x ∀ ( u , T ) ∈ H . (C 2) Proof.
Since T is a Hilbert space, the Riesz representation theorem guarantees thatthere exists a fixed temperature field θ ∈ T such that L{ T } = Ω ∇ θ · ∇ T d x . (C 3)Next, fix any pair ( u , T ) ∈ H and observe that, by virtue of the horizontal periodicity,the pair of translated fields ( u r,s , T r,s ) := ( u ( x + r, y + s, z ) , T ( x + r, y + s, z )) also belongsto H . By assumption S{ u , T } > L{ T } for all ( u , T ) ∈ H , so S{ u r,s , T r,s } > Ω ∇ θ · ∇ T r,s d x = Ω ∇ θ − r, − s · ∇ T d x , (C 4)where the second equality follows from a change of variables. The same change of variablesshows that S{ u r,s , T r,s } = S{ u , T } , so averaging (C 4) over all horizontal shifts r and s yields S{ u , T } > Ω ∇ θ · ∇ T d x = Ω θ ′ ( z ) ∂T∂z d x . (C 5)The expression on right-hand side of this inequality is a positive linear functional becauseit is the average of the positive functionals T → L{ T r,s } . Since the pair ( u , T ) ∈ H isarbitrary, inequality (C 2) follows upon setting q ( z ) := − θ ′ ( z ) and applying Lemma 2 toconclude that q is nondecreasing and square-integrable. Appendix D. Computational approach
The optimisation problems (4.12) and (4.25) can be discretised into SDPs followinga general strategy, and then solved using efficient algorithms for convex optimisation.This “discretise-then-optimise” approach preserves the linearity of the original infinite-dimensional problems and enables one to readily impose additional constraints, such asthe inequalities on ψ (0), ψ (1) and the monotonicity constraint on q , that are not easy toenforce following “optimise-then-discretise” strategies based on the numerical solution ofthe Euler–Lagrange equations for (4.12) and (4.25).The discretisation process starts by approximating the tunable functions ψ , q and theunknown fields ˆ T , ˆ T k and ˆ w k using a finite set of basis functions { Φ ( z ) , . . . , Φ n ( z ) } ,e.g., ψ ( z ) = n X i =1 A i Φ i ( z ) . (D 1)Here we use a single set of basis functions for simplicity, but different fields could beapproximated using different bases to improve accuracy or allow for varying degrees ofsmoothness. Note that while the functions ψ and q are arbitrary, so we are free to choose7such a finite-dimensional representation without much loss of generality, assuming thesame for the test functions ˆ T , ˆ T k and ˆ w k represents a relaxation of the constraintsin (4.12) and (4.25). Strictly speaking, therefore, our numerical results are not rigorousupper bounds on the vertical heat flux, but we expect convergence as n → ∞ .Substituting expansions such as (D 1) into the inequalities on S and ˜ S in (4.12)and (4.25) reduces them to quadratic polynomial inequalities, where the independentvariables are the (unknown) expansion coefficients of ˆ T , ˆ T k and ˆ w k , and the coefficientsdepend linearly on the optimisation variables—the scalars U , a , b and the expansioncoefficients of ψ and q . These quadratic polynomial inequalities are equivalent to positivesemidefinitess constraints on matrices that depend linearly on the polynomial coefficients,and hence on the optimisation variables. Moreover, the inequalities ψ (0) ψ (1) q ( z ) in (4.25) can be projected onto the expansionbasis to obtain a set of linear constraints on the expansions coefficients of ψ and q . Thediscrete problems are therefore SDPs (Boyd & Vandenberghe 2004) and can be solvedwith a variety of algorithms (see, e.g., Nemirovski 2006).To tackle (4.12) and (4.25), we used a finite-element approximation similar to that con-sidered by Fantuzzi et al. (2018). The reason for this choice is twofold. First, a piecewise-linear finite-element representation for q enables us to impose exactly the monotonicityconstraint, which is key to enforcing the minimum principle on the temperature asdiscussed in § ψ and q have steep boundary layers near the bottomboundary that cannot be approximated accurately at a reasonable computational costusing global polynomial expansions (e.g. Legendre series). Finite-element bases, instead,lead to SDPs with so-called chordal sparsity (Fukuda et al. sdpa-gmp (Fujisawa et al. et al. Appendix E. Estimates on the spectral constraint A R dependent variation of the parameters follows from enforcing the spectral con-straint. In the bulk of the domain the constraint can easily be satisfied for any profile ψ ( z ) such that ψ ′ ( z ) = a for z ∈ ( δ, − ε ). With this assumption, the only sign-indefiniteterm in (4.13) is ˆ ( a − ψ ′ ) ˆ w ˆ T d z = ˆ δ ( a − ψ ′ ) ˆ w ˆ T d z + ˆ − ε ( a − ψ ′ ) ˆ w ˆ T d z. (E 1)To estimate the integral over [0 , δ ], we apply the fundamental theorem of calculus, theboundary condition ˆ w k (0) = 0 and the Cauchy–Schwarz inequality to bound w ( z ) = ˆ z w ′ ( ξ ) d ξ √ z k w ′ k . (E 2)8Identical estimates show that T ( z ) √ z k T ′ k . Using these inequalities, we can estimate (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ˆ δ ( a − ψ ′ ) ˆ w ˆ T d z (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) k a − ψ ′ k L ∞ (0 ,δ ) ˆ δ | ˆ w | (cid:12)(cid:12) ˆ T (cid:12)(cid:12) d z δ k a − ψ ′ k L ∞ (0 ,δ ) k w ′ k k T ′ k . (E 3)Similar analysis near z = 1 yields (cid:12)(cid:12)(cid:12)(cid:12) ˆ − ε ( a − ψ ′ ) ˆ w ˆ T d z (cid:12)(cid:12)(cid:12)(cid:12) ε k a − ψ ′ k L ∞ (1 − ε, k w ′ k k T ′ k . (E 4)Substituting inequalities (E 3) and (E 4) into (4.13) shows that˜ S > aR k ˆ w ′ k − λ ( a, δ, ε, ψ ′ ) k ˆ w ′ k k ˆ T ′ k + b k ˆ T ′ k , (E 5)where λ ( a, δ, ε, ψ ′ ) := δ k a − ψ ′ k L ∞ (0 ,δ ) + ε k a − ψ ′ k L ∞ (1 − ε, . (E 6)The right-hand side of the last inequality is a homogeneous quadratic form in the variables k ˆ w ′ k and k ˆ T ′ k , so it is nonnegative if its discriminant is nonpositive, i.e., | λ ( a, δ, ε, ψ ′ ) | ab R . (E 7)Taking square roots and rearranging we conclude that ˜ S > δ k a − ψ ′ k L ∞ (0 ,δ ) + ε k a − ψ ′ k L ∞ (1 − ε, √ ab √ R . (E 8)This condition reduces to (5.3) when δ = ε and ψ is as in (5.1). Appendix F. Proof of Proposition 1
It is assumed that ( a, b, ψ, q ) ∈ A{ δ, ε } ∩ B{ γ } ∩ C{ δ, ε, R } . The first step is to bound ´ ψ ( z ) dz from above. To do this, we work back from z = 1, using the assumptions toestimate ψ .Let 1 − ε z
1. Since ψ ∈ C [0 ,
1] and ψ (1) = 0, the mean value theorem impliesthat there exist z ε ∈ (1 − ε,
1) such that − ψ ′ ( z c ) = ψ ( z )1 − z > ψ ( z ) ε . Since ( a, b, ψ ( z )) satisfy (E 8), it follows that ε k a − ψ ′ k L ∞ (1 − ε, a b R = ⇒ ε (cid:18) a + ψ ( z ) ε (cid:19) a b R = ⇒ ε (cid:16) ab (cid:17) + (cid:16) ab (cid:17) (cid:18) ε ψ ( z ) b − R (cid:19) + (cid:18) εψ ( z ) b (cid:19) ⇒ (cid:18) ε ψ ( z ) b − R (cid:19) > ε (cid:18) εψ ( z ) b (cid:19) = ⇒ ψ ( z ) ε b R (F 1)9Next, since ψ (1) = 0, the assumption that k ψ ′ k L ∞ (1 − ε, b gives | ψ ( z ) | ε k ψ ′ k L ∞ (1 − ε, bε , which in turn implies ε > b ψ ( z ) . (F 2)Since 1 − ε < z < ψ ′ ( z ) > z ∈ [ δ, − ε ] that k ψ k L ∞ ( δ, k ψ k L ∞ (1 − ε, bR − . (F 3)We now estimate ψ in the lower boundary layer [0 , δ ]. The mean value theorem impliesthat there exists z δ such that ψ ′ ( z δ ) = ψ ( δ ) − ψ (0) δ . Since ( a, b, ψ ( z )) satisfy (E 8), using the above equation and the assumption that 0 b ˆ ψ ( z ) dz δb ψ (0) + (1 − δ ) b k ψ k L ∞ ( δ, √ R (F 5)Next we consider the L component of the cost function. Using the assumption that ψ ′ , q are constant in an interval (1 / − γ, / γ ),14 b (cid:13)(cid:13) ψ ′ − q − b ( z − / (cid:13)(cid:13) > b ˆ
12 + γ − γ [ ψ ′ − q − b ( z − / d z > b ˆ + γ − γ ( z − / d z = bγ . (F 6)Combining the above estimate with (F 5) gives the stated result14 b (cid:13)(cid:13) ψ ′ − q − b ( z − / (cid:13)(cid:13) − ˆ ψ ( z ) dz > b γ − · √ R ! . REFERENCESBercovici, D.
Encyclopedia of Solid Earth Geophysics. Gupta,HK,(ed.) Springer . Bouillaut, Vincent, Lepot, Simon, Aumaˆıtre, S´ebastien & Gallet, Basile
Journalof Fluid Mechanics , R5.
Boyd, S. & Vandenberghe, L.
Convex Optimization . Cambridge University Press.
Chernyshenko, S. I. arXivpreprint arXiv:1704.02475 . Chernyshenko, S. I., Goulart, P. J., Huang, D. & Papachristodoulou, A.
PhilosophicalTransactions of the Royal Society A: Mathematical, Physical and Engineering Sciences (2020), 20130350.
Choffrut, A., Nobili, C. & Otto, F.
Journal of Differential Equations (4), 3860–3880.
Constantin, P. & Doering, C. R.
Physical Review E (4), 3192. Constantin, P. & Doering, C. R.
Nonlinearity (4), 1049–1060. Constantin, P. & Doering, C. R.
Journal ofStatistical Physics (1-2), 159–172. Debler, W. R.
Doering, C. R. & Constantin, P.
Physical Review E (5), 4087. Doering, C. R. & Constantin, P.
Physical Review E (6), 5957. Doering, C. R. & Constantin, P.
Journalof Fluid Mechanics , 263–296.
Doering, C. R. & Constantin, P.
Journal of Mathematical Physics (2), 784–795. Doering, C. R., Otto, F. & Reznikoff, M. G.
Journal of fluid mechanics ,229–241.
Emara, A. A. & Kulacki, F. A.
Journal of Heat Transfer (3), 531–537.
Fantuzzi, G., Goluskin, D., Huang, D. & Chernyshenko, S. I.
SIAMJournal on Applied Dynamical Systems (4), 1962–1988. Fantuzzi, G., Pershin, A. & Wynn, A.
Journal of Fluid Mechanics , 562–596.
Fantuzzi, G. & Wynn, A.
Physics Letters A (1-2), 23–32.
Fantuzzi, G. & Wynn, A.
Physical Review E (4), 043308. Fantuzzi, G., Wynn, A., Goulart, P. J. & Papachristodoulou, A.
IEEE Transactions onAutomatic Control (12), 6221–6236. Foias, C., Manley, O. & Temam, R.
Nonlinear Analysis: Theory, Methods &Applications (8), 939–967. Fujisawa, K., Fukuda, M., Kobayashi, K., Kojima, M., Nakata, K., Nakata, M. &Yamashita, M.
Tech. Rep. . Department of Mathematical and ComputingSciences, Tokyo Institute of Technology, Tokyo, Japan.
Fujisawa, Katsuki, Kim, Sunyoung, Kojima, Masakazu, Okamoto, Y. & Yamashita,Makoto
Tech. Rep. . Department of Mathematical andComputing Sciences, Tokyo Institute of Technology, Tokyo, Japan. Fukuda, M., Kojima, M., Murota, K. & Nakata, K.
SIAM Journalon Optimization (3), 647–674. Goluskin, D.
Journal of FluidMechanics , 36–56.
Goluskin, D.
Internally heated convection and Rayleigh-B´enard convection . Springer.
Goluskin, D. & Doering, C. R.
Journal of Fluid Mechanics , 370–386.
Goluskin, D. & Fantuzzi, G.
Nonlinearity (5), 1705–1730. Goluskin, D. & van der Poel, E. P.
Journal of Fluid Mechanics . Goluskin, D. & Spiegel, E. A.
Physics LettersA (1-2), 83–92.
Howard, L. N.
Journal of Fluid Mechanics (4),509–512. Jahn, M. & Reineke, H.-H.
Proceedings of the 5 th International Heat TransferConference , pp. 74–78. Tokyo.
Kakac, S., Aung, W. M. & Viskanta, R.
Washington, DC, Hemisphere Publishing Corp., 1985, 1191 p. . Kulacki, F. A. & Goldstein, R. J.
Journal of Fluid Mechanics (2), 271–287. Lee, S. D., Lee, J. K. & Suh, K. Y.
Journal of Heat Transfer (5), 679–682.
Lu, L., Doering, C. R. & Busse, F. H.
Journal of mathematical physics (7), 2967–2986. Malkus, W. V. R.
Proceedingsof the Royal Society of London. Series A. Mathematical and Physical Sciences (1161),196–212.
Mayinger, F., Jahn, M., Reineke, H. & Steibnberner, V.
Tech. Rep.
BMFT RS 48/1.Institut f¨ur Verfahrenstechnik der TU, Hanover Germany.
Miles, J. W.
Journal of Fluid Mechanics (4), 496–508. Nemirovski, A.
InternationalCongress of Mathematicians , , vol. 1, pp. 413–444.
Otero, J., Dontcheva, L. A., Johnston, H., Worthing, R. A., Kurganov, A., Petrova,G. & Doering, C. R.
Journal of Fluid Mechanics , 263–281.
Otero, J., Wittenberg, R. W., Worthing, R. A. & Doering, C. R.
Journal of Fluid Mechanics ,191–199.
Otto, F. & Seis, C.
Journal of mathematical physics (8), 083702. Peckover, R. S. & Hutchinson, I. H.
Physics of Fluids (7), 1369–1371. Plasting, S. C. & Kerswell, R. R.
Journal of Fluid Mechanics ,363–379.
Priestley, C. H. B.
Australian Journal of Physics (1), 202–209. Rosa, R. & Temam, R. M. arXiv preprint arXiv:2010.06730 . Spiegel, E. A.
TheAstrophysical Journal , 216.
Straus, J. M.
TheAstrophysical Journal , 179–189.
Tilgner, A.
Physical ReviewFluids (12), 123502. Tilgner, A.
Physical Review Fluids (1), 1–11. Tobasco, I., Goluskin, D. & Doering, C. R.
Physics Letters A (6), 382–386.
Tritton, D. J.
Nature (5522), 110–112.
Trowbridge, A. J., Melosh, H. J., Steckloff, J. K. & Freed, A. M.
Nature (7605), 79–81.
Tveitereid, M.
International Journal of Heat and Mass Transfer (3), 335–339. Waki, H., Nakata, M. & Muramatsu, M.
Computational Optimization and Applications (3), 823–844. Wang, Q., Lohse, D. & Shishkina, O.
Geophysical Research Letters , e2020GL091198. Wen, B., Chini, G. P., Dianati, N. & Doering, C. R.
PhysicsLetters A (41), 2931–2938.
Wen, B., Chini, G. P., Kerswell, R. R. & Doering, C. R.
Physical Review E (4), 043012. Whitehead, J. P. & Doering, C. R. a Internal heating driven convection at infinitePrandtl number.
Journal of mathematical physics (9), 093101. Whitehead, J. P. & Doering, C. R. b Ultimate State of Two-Dimensional Rayleigh–B´enard Convection between Free-Slip Fixed-Temperature Boundaries.
Physical ReviewLetters (24), 244501.
Whitehead, J. P. & Doering, C. R.
Journal of Fluid Mechanics , 241–259.
W¨orner, M., Schmidt, M. & Gr¨otzbach, G.
Journal of Hydraulic Research (6), 773–797. Yamashita, M., Fujisawa, K., Fukuda, M., Kobayashi, K., Nakata, K. & Nakata, M.
Handbook onsemidefinite, conic and polynomial optimization , pp. 687–713. Springer.
Yan, X.
Journal of Mathematical Physics45