Table Of Content

Principle of Maximum Entropy Applied to Rayleigh-Bénard Convection Takafumi Kita Department of Physics, Hokkaido University, Sapporo 060-0810, Japan (Dated: February 6, 2008) 7 A statistical-mechanical investigation is performed on Rayleigh-Bénard convection of a dilute 0 classicalgasstartingfromtheBoltzmannequation. Wefirstpresentamicroscopicderivationofbasic 0 hydrodynamicequationsandanexpressionofentropyappropriatefortheconvection. Thisincludes 2 an alternative justification for theOberbeck-Boussinesq approximation. We thencalculate entropy n change through the convective transition choosing mechanical quantities as independent variables. a AbovethecriticalRayleighnumber,thesystemisfoundtoevolvefromtheheat-conductinguniform J statetowardstheconvectiverollstatewithmonotonicincreaseofentropyontheaverage. Thus,the 6 principle of maximum entropy proposed for nonequilibrium steady states in a preceding paper [T. 1 Kita:J.Phys.Soc.Jpn.75(2006)114005]isindeedobeyedinthisprototypeexample. Theprinciple also providesa natural explanation for theenhancement of the Nusselt numberin convection. ] h c I. INTRODUCTION inaries, we perform a statistical-mechanical calculation e of entropy change through the convective transition by m In a preceding paper,1 we have proposed a principle choosingmechanicalquantitiesasindependentvariables. - t of maximum entropy for nonequilibrium steady states: It is worth pointing out that classical gases have been a The state which is realized most probably among possible used extensively for detailed experiments on Rayleigh- t s steadystateswithouttimeevolutionistheonethatmakes Bénard convection over the last two decades.5,8,9 Thus, . t entropy maximum as a function of mechanical variables. quantitativecomparisonsbetweentheoryandexperiment a are possible here. m We here apply it to Rayleigh-Bénardconvection of a dilute classical gas to confirm its validity. This paper is organized as follows. Section II derives - (i) equations of motion for the particle, momentum and d Rayleigh-Bénardconvectionisaprototypeofnonequi- n libriumsteadystates withpatternformation,andexten- energydensitiesand(ii)anexplicitexpressionforthedis- o sivestudieshavebeencarriedouttoclarifyit.2,3,4,5,6,7,8,9 tribution function f, both starting from the Boltzmann c equation. SectionIIIreducestheequationsof IIinaway However, no calculation seems to have been performed [ § on entropy change through the nonequilibrium “phase appropriateto treat Rayleigh-Bénardconvectionof a di- 3 transition,” despite the fact that entropy is the key luteclassicalgas. Thisincludesasystematicjustification v oftheOberbeck-Boussinesqapproximationandaderiva- concept in equilibrium thermodynamics and statistical 1 tion of the expression of entropy for convection. Section mechanics. There may be at least four reasons for 7 IVtransformstheequationsof IIIintothosesuitablefor it. First, there seems to have been no established ex- 2 § periodic structures with the stress-free boundaries. Sec- 1 pression of nonequilibrium entropy. Second, the stan- tionVpresentsnumericalresultsobtainedbysolvingthe 1 dard starting point to describe Rayleigh-Bénardconvec- 6 tion is a set of deterministic equations for the particle, equationsof IV. Itisshownexplicitlythattheprinciple § 0 momentum and energy flows, with all the thermody- ofmaximumentropyisindeedobeyedbythe convection. / Concluding remarks are given in VI. t namic effects pushed into phenomenological parameters § a of the equations.2,6 Third, one additionally adopts the m Oberbeck-Boussinesq approximation to the equations in - the conventional treatment.2,6,10 Despite many theoreti- II. DISTRIBUTION FUNCTION, d n cal efforts overa long period,11,12,13,14,15 a well-accepted CONSERVATION LAWS AND ENTROPY o systematic justification for it seems still absent, thereby c preventing a quantitative estimation of entropy change. A. The Boltzmann equation and entropy v: Fourth,thereisambiguityonwhattochooseasindepen- i dent variables of entropy for open systems. We shallconsidera monatomic dilute classicalgasun- X In a preceding paper,1 we have derived an expression der gravity. This system may be described by the Boltz- ar of nonequilibrium entropy together with the evolution mann equation:16,17 equations for interacting bosons/fermions. We here ap- plythemtoaclassicalgasofthedilutehigh-temperature ∂f p ∂f ∂f + mg = . (1) limit where the evolution equations reduce to the Boltz- ∂t m · ∂r − ∂p C z mannequation. Wecarryoutamicroscopicderivationof the hydrodynamicequationsforthe particle,momentum Heref=f(p,r,t)isthedistributionfunction,tthetime, and energy densities (i.e., the basic conservation laws) p the momentum, m the mass, r the space coordinate, fromtheBoltzmannequation. Wethenprovideasystem- and g the acceleration of gravity. With a unified de- atic justification for the Oberbeck-Boussinesq approxi- scription of classical and quantum statistical mechanics mation to describe the convection. With these prelim- in mind, we choose the normalizationoff such thatit is 2 dimensionless and varies between 0 and 1 for fermions. B. Conservation laws The collision integral is given explicitly by C We next consider conservation laws which originate (p,r,t) C from the Boltzmann equation. The basic physical quan- =~2 (2dπ3p~1)3 (2dπ3~p)′3 (2dπ3p~′1)3|Vp′−p|2(fp′fp′1−fpfp1) tvi(trie,st)r,eltehveatnetmtopetrhaetmuraerTe(trh,etd),enthsietymno(mre,nt)t,utmheflvuexlodceinty- Z Z Z ×(2π)4δ(Ep′+Ep′1−Ep−Ep1)δ(p′+p′1−p−p1), (2) lsoitcyaltevnesloocriΠty(rv,,ta)nindtthheerheefaerteflnucxefdreanmseitymjovi(nrg,tw).itThhtheye Q whereV isFouriertransformoftheinteractionpotential are defined by q and E =p2/2m. p d3p Equation (1) also results from eq. (63) of ref. 1 as fol- n(r,t)= f(p,r,t), (7a) lows: (i)approximatethespectralfunctionasA(pε,rt)= (2π~)3 Z 2πδ(ε E mgz); (ii) substitute the second-order self- p − − energy of eq. (66) into the collision integral of eq. (64); 1 d3p p (iii) take the high-temperature limit; and (iv) integrate v(r,t)= f(p,r,t), (7b) eq. (63) over ε to obtain an equation for n(r,t) (2π~)3m Z ∞ dε f(p,r,t) A(pε,rt)φ(pε,rt). (3) 2 d3p p¯2 ≡Z−∞2π T(r,t)= 3k n(r,t) (2π~)32mf(p,r,t), (7c) B Z Thiswholeprocedureamountstotreatingtheinteraction potential only as the source of dissipation in the dilute high-temperaturelimit,thusneglectingcompletelyitsin- d3p p¯p¯ Π(r,t)= f(p,r,t), (7d) fluence on the density of states, i.e., on the real part of (2π~)3 m Z the self-energy. With a change of variables p =p+q, p′=P q′/2 1 and p′ = P +q′/2 in eq. (2), the collision integ−ral is d3p p¯2 p¯ 1 j (r,t)= f(p,r,t), (7e) transformed into Q (2π~)32mm Z d3q q C = (2π~)3 m dσ[fp+(q−q′)/2fp+(q+q′)/2−fpfp+q], fwroitmh pe¯q≡s.p(C−1ma)v,.(CT1hbe)e,x(pCre2s0saio),ns(C(71a4))-(a7ned) a(Clso20rbes)uoltf Z Z (4) ref. 1, respec·tively, b·y notin·g eq. (3)·, neglecting·the in- where dσ=dΩq′ dq′δ(q′ q) mV(q′−q)/2 /4π~2 2 is the teractionterms,andidentifying thelocalinternal-energy differentialcrosssectioni−nthe ce|nter-of-m|asscoordinate density ãs ˜=3nk T. R (cid:2) (cid:3) E E 2 B of the scattering with dΩq′ denoting the infinitesimal To obtain the number, momentum and energy conser- solid angle. We shall use the contact interaction with vation laws, let us multiply eq. (1) by 1, p and p¯2/2m, no q dependence in V , where dσ reduces to respectively, and perform integration over p. The con- q tribution from the collision integral (2) vanishes in all dσ =a2dΩq′ dq′δ(q′−q), (5) tehrgeythcorenesecravsaetsiodnusethtroouthgehpthaertcicollel,ismionom.16e,n17tuTmheanredsuelnt-- Z ing hydrodynamic equations can be written in terms of witha mV /4π~2denotingthescatteringlength. Now, the quantities of eq. (7) as ≡ | | thedifferentialcrosssectionacquirestheformofthetwo- particle collision between hard-sphere particles with ra- ∂n dius a.18 +∇(nv)=0, (8a) ∂t Entropy per unit volume is given in terms of f by k d3rd3p ∂v 1 S = B f(logf 1), (6) +v ∇v+ ∇Π+ge =0, (8b) − (2π~)3 − ∂t · mn z V Z with denotingthevolume. Thisexpressionalsoresults fromVeq.(69)ofref.1forentropydensityby: (i)adopting 3nk ∂T +v ∇T +∇ j +Π:∇v =0, (8c) B Q the quasiparticle approximationA(pε,rt)=2πδ(ε E 2 ∂t · · − p− (cid:18) (cid:19) mgz); (ii) defining f by eq. (3) above; (iii) taking the high-temperature limit; and (iv) performing integration where ez is the unit vector along the z axis and A:B ionvetrher.BTolhtzemtearnmnH−1-fuinnctthieonin16tebgurtanndatoufraeqll.y(r6e)suisltasbfsreonmt edqeunaottieosntshaereteindseonrticparloidnucfotr:mAw:iBth≡eqs.i(jCAi2j)B,j(iC. 9T)haensde P · · the procedure (iii) above. Indeed, eq. (6) reproduces the (C21) of ref. 1, respectively, with U = mgz and ˜ = correct expression of entropy in equilibrium.19 3n·k T. E 2 B 3 C. The Enskog series with P=nk T the pressure and 1 the unit tensor. The B left-hand side of eq. (12) is thereby transformed into We nowreducethe wholeproceduresofsolvingeq.(1) ∂f(eq) p ∂f(eq) p¯ g for(p,r,t)tothoseofsolvingeq.(8)for(r,t). Weadopt + + z f(eq) thewell-knownEnskogmethod16,17 forthispurpose,i.e., ∂t m · ∂r kBT the expansion from the local equilibrium. We here de- k¯2 5 p¯ =f(eq) 2 k¯k¯ 1 :∇v+ k¯2 ∇lnT , (15) scribe the transformation to the extent necessary for a − 3 −2 m· (cid:20) (cid:18) (cid:19) (cid:18) (cid:19) (cid:21) later application to Rayleigh-Bénardconvection. Let us expand the distribution function formally as with k¯ a dimensionless quantity defined by f(p,r,t)=f(eq)(p,r,t) 1+ϕ(1)(p,r,t)+ , (9) k¯ p¯/ 2mkBT . (16a) ≡ ··· h i Itisconvenienttointroducpetheadditionaldimensionless wheref(eq) isthelocal-equilibriumdistributiongivenex- quantities: plicitly by m vˆ v, rˆ n1/3r, (16b) (2π~)3n p¯2 ≡ 2k T ≡ f(eq) = exp , (10) r B (2πmk T)3/2 −2mk T B (cid:18) B (cid:19) and the mean-free path: with p¯ p mv. This f(eq) has been chosen so as 1 to satisf≡y the−two conditions:16,17 (i) the local equilib- l . (17) ≡ 4√2πa2n rium condition that the collision integral vanishes; (ii) eqs. (7a)-(7c) by itself. It hence follows that the higher- Using eqs. (16) and (17) and noting eqs. (5), (10), (13) ordercorrectionsϕ(j) (j=1,2, )ineq.(9)shouldobey and (15), we can transform eq. (12) into the dimension- ··· the constraints: less form: d3p p¯nf(eq)ϕ(j) =0 (n=0,1,2). (11) 2√2 d3q d3q′δ(q′−q)e−k¯2−(k¯+q)2 (2π~)3 √πln1/3 4π 4π q Z Z Z ϕ(1) +ϕ(1) ϕ(1) ϕ(1) Note that we have incorporated five space-time depen- × k+(q−q′)/2 k+(q+q′)/2− k − k+q dbeendteptearrmamineetdercsominplfe(teeql)y, wi.eit.h, neq,sv. (a8na)d-(T8c,)wohficcohnscearn- =e−hk¯2 2 k¯k¯ k¯21 :∇ˆvˆ+ k¯2 5 k¯ ∇ˆ lniT , (18) − 3 −2 · vation laws. The remaining task here is to express the (cid:20) (cid:18) (cid:19) (cid:18) (cid:19) (cid:21) extra quantities Π and jQ in eq. (8) as functionals of n, where∇ˆ ∂/∂rˆ,andwehaveredefinedϕ(1) asafunction v and T. ≡ of k p/√2mk T. Similarly, eq. (11) now reads B Letussubstitute eq.(9)intoeq.(1), regardthespace- ≡ time differential operators on the left-hand side as first- d3ke−k¯2k¯nϕ(1) =0 (n=0,1,2). (19) order quantities, and make use of the fact that the colli- k sion integral (2) vanishes for f(eq). We thereby arrive at Z the first-order equation: Equation(18)withsubsidiarycondition(19)formsalin- ear integral equation for ϕ(1). k ∂f(eq) + p ∂f(eq) + p¯zg f(eq) = (1), (12) The right-hand side of eq. (18) suggests that we may ∂t m · ∂r k T C seek the solution in the form:16,17 B k¯2 where (1) is obtained from eq. (4) as ϕ(1) = ln1/3 A5/2(k¯2) k¯k¯ 1 :∇ˆvˆ C k − − 3 (cid:20) (cid:18) (cid:19) (1) = d3q q dσf(eq)f(eq) +A3/2(k¯2)k¯ ∇ˆ lnT , (20) C (2π~)3 m p p+q · Z Z i ϕ(1) +ϕ(1) ϕ(1) ϕ(1) . (13) whereA5/2 andA3/2 aretwounknownfunctions;theuse × p+(q−q′)/2 p+(q+q′)/2− p − p+q of fractions α=5/2 and 3/2 to distinguish them will be h i rationalizedshortly. Substituting eq.(20)into it, we can With eq. (10), the derivatives of f(eq) in eq. (12) are transform eq. (18) into separate equations for Aα as transformedintothoseofn,vandT. Wethenremovethe taipmperodxeirmivaattiiovne,swbhyeuresing eq. (8) in the local-equilibrium 2√√π2 d43πq d43πq′ δ(q′q−q)e−k¯2−(k¯+q)2 Tk¯α+Tk¯α+q Z Z h Π(eq)=P1, jQ(eq)=0, (14) −Tk¯α+(q−q′)/2−Tk¯α+(q+q′)/2 =e−k¯2Rαk¯, (21) i 4 where tensor α and α are defined in Table I together The quantities Iα (k,q,q) andIα (k,q, q) areobtained withanothertRenksor Tαk. Sincethe factore−k¯2 ispresent fromeq.(27)byℓrℓe′movingthe inℓtℓe′gralo−verdΩq′/4π and on the right-hand sidWekof eq. (21), we expand Aα(ε) fur- settingq′ qandq′ qintheintegrand,respectively. → →− ther in the Sonine polynomials Sα(ε) as16,17 The first few series of eq. (26) are easily calculated ℓ analytically as 5/2 = 3/2 = 1, 5/2 = 3/2 = 1/4, ∞ T00 T11 T01 T12 − Aα(ε)= cαℓSℓα(ε). (22) T151/2 = 205/48 and T232/2 = 45/16, in agreement with the values given below eq. (10.21,3) of Chapman and ℓ=0 X Cowling.16 Thematrixelementforageneralℓℓ′ caneval- Uonsleyoffotrwcoondviffenerieennctecotomtprlaentesfsoertms {tShℓαe}rℓig(hαt=-h5a/n2d,3si/d2e)oisf u(2a4t)edisnsoulmveedricbayllcyu.ttWinigththTeℓαℓi′nfiannditeRsαℓeritehsuastoabfitaniinteedv,aleuqe. eq.(21)intoavectorwithasinglenonzeroelement. With ℓ , and ℓ is increased subsequently to check the conver- c c eq.(20)and(22)andthe orthogonalityofSα(ε), wefind gence. Table II lists the values of cα thereby obtained. ℓ ℓ that the constraint (19) reduces to the single condition: Thoseofc5/2 andc3/2 are about2%largerinmagnitude 0 1 thantheanalyticones5√π/4and 15√π/16withℓ =0 c3/2 =0. (23) − c 0 and 1 for α=5/2 and 3/2, respectively. This rapid con- vergence as a function of ℓ was already pointed out by We hence remove the ℓ=0 term of α=3/2 from eq. (22) c Chapman and Cowling.16 in the subsequent discussion. We now take the tensor Substituting eqs. (9), (10) and (20) into eqs. (7d) and (α=5/2)orthevector(α=3/2)productofeq.(21)with Sα(k¯2) α/4π andperformintegrationoverk¯. Equation (7e), wearriveatthefirst-ordercontributionstothe mo- ℓ Wk¯ mentum flux density tensorandthe thermalflux density (21)istherebytransformedintoanalgebraicequationfor the expansion coefficients cα as as { ℓ}ℓ ∂v ∂v 2 ℓ′ Tℓαℓ′cαℓ′ =Rαℓ . (24) Π(ij1) =−mnν(cid:18)∂rji+∂rji −δij3∇·v(cid:19), (28a) X Here α is defined by Rℓ j(1) = 3nk κ∂T , (28b) α d3ke−k2Sα(k2)( α, α) Q −2 B ∂r Rℓ ≡ 4π ℓ Wk Rk respectively, where ν and κ are the kinematic viscosity Z 5 and the thermal diffusivity,2,6 respectively, defined by √πδ :α=5/2 ℓ0 4 = , (25)  −1156√πδℓ1 :α=3/2 ν = 14lr2kmBTc50/2, (29a) with ( , )  : and for α= 5/2 and 3/2, raebslpeek¯ctWivekRly.qA≡/l2sWoin, TtℓRhαℓ′eik¯s oinbttWeaginr·aReldawsith a changeofvari- κ=−56l 2kmBTc31/2. (29b) → − r α Thesequantitiesclearlyhavethe samedimension. Using Tℓℓ′ = 2√2 ∞dke−2k2k2 ∞dqe−q2/2q3 ∞dq′δ(q′ q) tahnedmcoansswtaenlltavsotluhmesepCeci,ficwehecaatnaitntcroondsutacentapnriemsspuorretaCnPt V √π − Z0 Z0 Z0 dimensionless quantity Pr = (ν/κ)(CP/CV) called the × Iℓαℓ′(k,q,q)+Iℓαℓ′(k,q,−q)−2Iℓαℓ′(k,q,q′) , (26) Prandtl number.16 Adopting CP/CV =5/3 of the ideal monatomicgas,wefindPr=0.66fromeq.(29)andTable with(cid:2)Iα (k,q,q′) defined by (cid:3) ℓℓ′ II, which is in excellent agreement with the value 0.67 Iℓαℓ′(k,q,q′) ftohrerAmrodatynTam=ic27a3nKd.1t7raTnasbploertIIcIoleiffistcsievnatlsueasrooufnrdelervoaonmt ≡ d4Ωπq d4Ωπq′Sℓα((k−q/2)2)Sℓα′((k−q′/2)2) temTpheursa,twuerehaavte1sautcmcesfosfruNllye,exAprreasnsdedaΠir.andj interms Z Z Q ×(Wkα−q/2,Wkα−q′/2). (27) ofn, v andT as eqs.(14), (28) and(29) within the first- order gradient expansion. Now, eq. (8) with eqs. (14), TABLE I: Quantities appearing in eq. (21). Here Sα(ε)=1 0 and S1α(ε)=1+α−ε. TABLE II:Values of cαℓ obtained by solving eq. (24). α Wkα Rαk Tkα α cα0 cα1 cα2 cα3 cα4 5/2 kk (k2/3)1 2Sα(k2) α Aα(k2) α 5/2 2.2511 0.1390 0.0233 0.0058 0.0018 − 0 Wk Wk 3/2 k Sα(k2) α Aα(k2) α 3/2 0 1.7036 0.1626 0.0371 0.0117 − 1 Wk Wk − − − − 5 (28) and (29) forms a closed set of equations for the five the entropy change through the convective transition on parameters n, v and T incorporated in feq(p,r,t). Af- a firm ground. tersolvingthem,wecanobtainthe distributionfunction f(p,r,t) by eqs. (9), (10), (20), (22) and Table II, and subsequently calculate entropy by eq. (6). A. Introduction of dimensionless units We first introduce a characteristic temperature T¯ de- III. APPLICATION TO RAYLEIGH-BEŃARD fined by CONVECTION T¯ 2¯/3n¯k . (31) B ≡ E We now apply the equations of II for n, v and We then adopt the units where the length, velocity and T to Rayleigh-Bénard convection of§a dilute classical energy are measured by d, k T¯/m and k T¯, respec- B B monatomic gas confined in the region d/2 z d/2 tively. Accordingly,wecarryoutachangeofvariablesas − ≤ ≤ p and L/2 x,y L/2. The gas is heated from below so − ≤ ≤ that m t=d t′, r =dr′, (32a) k T¯ T(x,y,z= d/2)=T0 ∆T/2, ∆T >0. (30) r B ± ∓ and The thickness d and the lateral width L are chosen as l d L. It hence follows that (i) there are enough n′ k T¯ ≪ ≪ n= , v = B v′, T =T¯T′. (32b) collisions along z and (ii) any effects from the side walls d3 m r may be neglected. We eventually impose the periodic Letussubstituteeqs.(14),(28)and(29)intoeq.(8)and boundary condition in the xy plane. subsequently perform the above change of variables. We We study this systemby fixing the totalparticle num- thereby obtain the dimensionless conservation laws: ber, total energy, and total heat flux through z= d/2. − This is equivalent to choosing the average particle den- ∂n′ sity n¯, the average energy density ¯, and the average +∇′ (n′v′)=0, (33a) heat flux density ¯j at z = d/2 aEs independent vari- ∂t′ · Q ables; hence T =T (n¯, ¯,¯j )−and ∆T =∆T(n¯, ¯,¯j ) in 0 0 E Q E Q ∂v′ eq. (30). The latter two conditions also imply, due to n′ +n′v′ ∇′v′+∇′P′ ′[n′ν′(∇′v′+ ′v′)] theenergyconservationlaw,thatthereisaverageenergy ∂t′ · − ∇i i ∇i i tflouxz.dTenhseitfyac¯jtQjutshtirfioeusgohuarncyhocircoesosfs¯jeQctaiosnanpeinrpdeenpdenicduelnatr +32∇′(n′ν′∇′·v′)+n′Ug′ez =X0, (33b) variable to specify the system. It should be noted that this energy flow in the container may be due partly to a ∂T′ 1 2 macroscopic motion of the gas. +v′ ∇′T′ ∇′ (n′κ′∇′T′)+ T′∇′ v′ ∂t′ · − n′ · 3 · The standard theoretical treatment of Rayleigh- 2 BénardconvectionstartsfromintroducingtheOberbeck- 2ν′ 1 ∂v′ ∂v′ 2 TBo.2uHssoinweesvqera,ptphriosxaimppartoioxnimtaotitohneseeqeumastinoontstfoorhnav,evbaenend − 3  ij 2 ∂rj′i + ∂rji′! − 3(∇′·v′)2=0, (33c) X justified in a widely accepted way. We here develop a   with P′=n′T′ and systematicapproximationschemefortheequationsin II § appropriate to treat Rayleigh-Bénard convection, which ν m κ m mgd ν′ , κ′ , U′ . (34) will be shown to yield the equations with the Oberbeck- ≡ d k T¯ ≡ d k T¯ g ≡ k T¯ Boussinesq approximation as the lowest-order approxi- r B r B B mation. This consideration also enables us to estimate Animportantdimensionlessquantityofthesystemisthe Rayleigh number R defined by U′∆T′ gT¯−1∆Td3 g TABLEIII:Valuesofrelevantthermodynamicandtransport R = , (35) ≡ ν′κ′ νκ coefficients under atmospheric pressure given in CGS units [α V−1(∂V/∂T)]. where T¯−1 appears as the thermal expansion coefficient ≡ α of the ideal gas. ◦ ◦ ◦ Ne(0 C) Ar(0 C) air (0 C) The above equations will be solved by fixing n¯, ¯= mαn 03..96060×1100−−33 13..7687×1100−−33 13..2697×1100−−33 3n¯kBT¯/2 and ¯jQ, as already mentioned. These conEdi- ν 33.0×10−2 11.8×10−2 13.2×10−2 tions are expressed in the dimensionless form as κ 83.1×10−2 29.2×10−2 25.9×10−2 1 × × × n′(r′)d3r′ =n¯′, (36a) L′2 Z 6 1 3P′+ 1n′v′2+n′U′z′ d3r′ = 3n¯′, (36b) we first expand n and T as L′2 2 2 g 2 Z (cid:18) (cid:19) ∞ ∞ n=n¯ 1+ nˆ(ℓ) , T =1+ T(ℓ). (40a) 1 3 ∂T′ ℓ=1 ! ℓ=1 dx′ dy′ n′κ′ =¯j′ , (36c) X X −L′2 Z Z 2 ∂z′(cid:12)z′=−12 Q Witheqs.(29),(34)and(38),wenextexpanddimension- (cid:12)(cid:12) less parameters ν =(l/4)√2Tc5/2, κ= (5l/6)√2Tc3/2 where P′=n′T′, and integrations(cid:12) extend over L′/2 0 − 1 x′,y′ L′/2 and 1/2 z′ 1/2. Equation (3−6b) ha≤s and Ug as been ≤obtained by−integr≤ation≤of (p2/2m+mgz)f over r ∞ ∞ and p with eq. (7), whereas eq. (36c) originatesfrom eq. ν = ν(ℓ), κ= κ(ℓ), Ug =Ug(2), (40b) (28b). ℓ=2 ℓ=2 X X To make an order-of-magnitude estimate for the parameters in eqs. (34) and (35), consider Ar of 273K at where ν(2) =(l/4)√2c50/2 and κ(2) =−(5l/6)√2c31/2 are 1atm confined in a horizontal space of dcm with the constants with l=l(n¯). It also follows from ∆T δ and ∼ temperature difference ∆T K. Using Table II, we then eq. (36c) that obtain the numbers: ¯j =¯j(3). (40c) Q Q ν′ =4.95 10−6/d, κ′ =1.22 10−5/d, × × (37a) It remains to attach ordersof magnitude to the differen- U′ =1.73 10−6d, ∆T′ =3.67 10−3∆T, g × × tial operators and j=nv. In this context, we notice that the Oberbeck-Boussinesq approximationyields a critical and RayleighnumberR whichisingoodquantitativeagree- c R=1.04 102d3∆T. (37b) ment with experiment.7 The fact tells us that the pro- × cedure to attach the orders should be carried out so as to reproduceR of the Oberbeck-Boussinesqapproxima- The criticalRayleighnumber R forthe convectivetran- c c sition is of the order 103,2 which is realized for d 2cm tion. The requirement yields ∼ and∆T 1K.Wenowobservethatthedimensionlesspa- ∞ ∼ ∂ rametershavethefollowingordersofmagnitudeinterms j =n¯ jˆ(ℓ+1/2), =O(δ1.5), ∇=O(δ−0.25). of δ 10−3: ∂t ≡ Xℓ=1 (40d) ν′,κ′,Ug′ ∼δ2, ∆T′ ∼δ, R∼δ−1. (38) See eq. (58) below and the subsequent comments for details. The above power-counting scheme will be shown Thus, Rayleigh-Bénard convection is a phenomenon to provide not only a justification of the Oberbeck- where two orders of magnitude (i.e., δ and δ2) are rele- Boussinesq approximation but also a systematic treat- vant. Fromnowonweshalldropprimesineveryquantity ment to go beyond it. of eqs. (33) and (36). Letussubstitute eq.(40)toeqs.(33b)and(36b). The contributions of O(δ) in these equations read P(1)=0 ∇ and P(1)d3r=0,respectively,withP(1)=n¯(nˆ(1)+T(1)). B. Omission of the number conservation law We hence conclude P(1)=0, i.e., R Let us write eq. (33) in terms of j=nv instead of v. nˆ(1) = T(1). (41) − It follows from the vector analysis that the vector field j can be written generally as j=∇Φ+∇ A, where Φ It also follows from eqs. (36a) and (36c) with eqs. (40) and A correspondsto the scalarandvector×potentials of and (41) that T(1) should obey theelectromagneticfields,respectively. Wethenfocusin thefollowingonlyonthosephenomenawherethecurrent T(1)d3r =0, (42a) density satisfies ∇Φ=0, i.e., Z ∇·j=0. (39) 1 L/2 L/2 ∂T(1) 2¯j(3) dx dy = Q . (42b) This implies thatwe maydropeq.(33a)fromeq.(33) to L2 Z−L/2 Z−L/2 ∂z (cid:12)z=−1/2 −3n¯κ(2) treat only eqs. (33b) and (33c). (cid:12) (cid:12) Equation (42) is still not suffi(cid:12)cient to determine T(1). It turns out below that the required equation results from C. Expansion in δ the O(δ3) and O(δ2.5) contributions of eqs. (33b) and (33c), respectively. Equation (38) suggests that we may solve eqs. (33b) Next,collectingtermsofO(δ2)ineqs.(33b)and(36b) and (33c) in powers of δ. Noting eqs. (36a) and (36b), yield P(2)+n¯U(2)e =0 and (3P(2)+n¯U(2)z)d3r = ∇ g z 2 g R 7 0, respectively. Hence P(2) = n¯U(2)z. Noting P(2) = D. Expression of entropy g n¯(nˆ(2)+nˆ(1)T(1)+T(2)) and usi−ng eq. (41), we obtain We now write down the expression of entropy in pow- nˆ(2) =(T(1))2 T(2) U(2)z. (43) − − g ers of δ. Entropy of the system can be calculated by It follows from eqs. (36a) and (36c) with eqs. (40) and eq. (6), where the distribution function f is given by (43) that T(2) should obey eq. (9) with eqs. (10) and (20). Hereafter we shall drop the superscript in ϕ(1), which specifies the order in the [(T(1))2 T(2)]d3r=0, (44a) gradient expansion, to remove possible confusion with − Z the expansion of eq. (40). Thus, f is now expressed as f=f(eq)(1+ϕ). L/2 L/2 ∂T(2) T(1)∂T(1) We first focus on ϕ and write eq. (20) in the present dx dy + =0. Z−L/2 Z−L/2 (cid:18) ∂z 2 ∂z (cid:19)(cid:12)z=−1/2 unitswithnotingeq.(16). We thenrealizethatϕispro- (cid:12) (44b) portional to l∇T or l∇v, which are quantities of O(δ3) (cid:12) In deriving eq. (44b), use has been made(cid:12)of (nκ)(3) = and O(δ3.5) in the expansion scheme of eq. (40), respec- n¯κ(2)T(1)/2whichresultsfromκ lT1/2 andl n−1; see tively. It hence follows that there is no contribution of eqs. (17) and (29b). Equation (4∝4) forms cons∝traints on O(δ2) from ϕ. In contrast, f(eq) yields terms of O(δ2), the higher-order contribution T(2), which will be irrele- as seen below. Thus, we only need to consider f(eq). vant in the present study, however. Let us write f(eq) of eq.(10) in the presentunits, sub- Finally,wecollecttermsofO(δ3)ineq.(33b)toobtain stitute eq. (40) into it, and expand the resulting expression in powers of δ. We also drop terms connected with ∂jˆ(1.5) +jˆ(1.5) ∇jˆ(1.5)+ ∇P(3) j (i.e., v) which have vanishing contribution to S within ∂t · n¯ O(δ2) after the momentum integration in eq. (6). We ν(2) 2jˆ(1.5)+nˆ(1)U(2)e =0, (45) thereby obtain the relevant expansion: − ∇ g z wherewehaveusedeq.(39). Wefurtheroperate∇×∇× (2π~)3n¯ 2 totheaboveequationandsubstituteeq.(41). Thisyields f(eq) = e−ε 1+ nˆ(ℓ) (2π)3/2 ! ℓ=1 X ∂ 2 −∂t∇2jˆ(1.5)+∇×∇×(jˆ(1.5)·∇jˆ(1.5)) 1+u(1)(ε) T(ℓ)+u(2)(ε)(T(1))2 , (47) ×" # +ν(2)(∇2)2jˆ(1.5)+Ug(2)(ez∇2−ez·∇∇)T(1) =0.(46a) Xℓ=1 On the other hand, terms of O(δ2.5) in eq. (33c) lead to where ε=p2/2, and u(1) and u(2) are defined by ∂T(1) 3 1 15 +jˆ(1.5) ∇T(1) κ(2) 2T(1) =0. (46b) u(1)(ε)=ε , u(2)(ε)= ε2 5ε+ . (48) ∂t · − ∇ − 2 2 − 4 (cid:18) (cid:19) Equation (46) forms a set of coupled differential equa- Let us substitute eq. (47) into eq. (6) and carry out in- tions for T(1) and jˆ(1.5), which shouldbe solvedwith eq. tegration over p. The contribution of O(1) is easily ob- (42). Itisalmostidenticalinformwiththatderivedwith tained as (k =1) the Oberbeck-Boussinesq approximation, predicting the B samecriticalRayleighnumberR aswillbeshownbelow. c (2π)3/2 5 ThewholeconsiderationsonRayleigh-Bénardconvection S(0)=n¯ ln + , (49) (2π~)3n¯ 2 presented in the following will be based on eq. (46) with (cid:20) (cid:21) eq. (42). which is just the equilibrium expression19 for density n¯ Two comments are in order before closing the subsec- andtemperatureT¯intheconventionalunits,asitshould. tion. First, if we apply the procedure of deriving eq. Next, we find S(1)=0 due to eqs. (41) and (42a). Thus, (46) to the O(δ4) andO(δ3.5) contributionsof eqs.(33b) the contribution characteristic of heat conduction starts and (33c), respectively, we obtain coupled equations for from the second order. A straightforward calculation the next-order quantities T(2) and jˆ(2.5), which should yields be solved with eq. (44). Thus, we can treat higher- osirodnerscchoenmtreib.uStieocnosnsdy,seteqm. (a4t5ic)almlyayinbteherepgraersdenedt eaxsptahne- S(2) = 5 T(1) 2d3r. (50) −4 equation to determine P(3) for given T(1) and jˆ(1.5). It V Z (cid:0) (cid:1) yields a relation between T(3) and nˆ(3), which in turn Equation(50) is the basic starting point to calculate the leads to the constraint for T(3) upon substitution into entropy change through the convective transition. Note eqs. (36a) and (36c). On the other hand, the equation that we have fixed ¯j in the present consideration, i.e., Q for T(3) originates from the O(δ5) and O(δ4.5) contribu- the initial temperature slope as seen from eq. (42b). tionsofeqs.(33b)and(33c),respectively. Now,onemay Atthisstate,itmaybeworthwhiletopresentaqualita- understand the hierarchy of the approximation clearly. tiveargumentonentropyofRayleigh-Bénardconvection. 8 With the initialtemperature slope fixedas eq.(42b), eq. B. Instability of the conducting state (50) tells us that entropy will be larger as the temperature profile becomes more uniform between z= 1/2. We next check stability of the conducting solution by ± The conducting state with v=0 has the linear tempera- adding a small perturbation given by tureprofile,asshownshortlybelowineq.(52). Thus,any increase of entropy over this conducting state is brought jˆ(1.5)(r,t)=eλteik⊥·r(δjssink ζ+δjc cosk ζ), (54a) z ⊥ z about by reducing the temperature difference between z= 1/2. Such a state necessarily accompanies a tempera±ture variation which is weaker around z = 0 than T(1)(r,t)= ∆Thcz+δTeλteik⊥·rsinkzζ, (54b) − near the boundaries z= 1/2. This temperature profile ± where ζ z+1/2, k = ℓ π (ℓ = 1,2, ) from eq. isindeedanessentialfeature ofRayleigh-Bénardconvec- (51), and≡δjc denotesz a ve3ctor 3in the xy···plane. Let tionwhichshowsupasanincreaseoftheNusseltnumber ⊥ us substitute eq. (54) into eq. (46) and linearize it with (i.e.,theefficiencyofheattransport)underfixedtemper- respect to the perturbation. This leads to ature difference.2,7 Combining eqs. (42b) and (50) with the experimentalobservationon the Nusselt number, we (λ+ν(2)k2)k2δjs U(2)k2δTe =0, (55a) thereby conclude without any detailed calculations that − g ⊥ z entropy of Rayleigh-Bénard convection should be larger thanentropyoftheconductingstateinthepresentcondi- (λ+ν(2)k2)k2δjc iU(2)k k δT =0, (55b) tionswith¯j =const. Thus,Rayleigh-Bénardconvection ⊥− g ⊥ z Q is expected to satisfy the principle of maximum entropy given at the beginning of the paper. We shall confirm (λ+κ(2)k2)δT ∆T δjs =0. (55c) − hc z this fact below through detailed numerical studies. The components δjs and δjc in the xy plane are ob- ⊥ ⊥ tained from eqs. (55a) and (55b) as IV. PERIODIC SOLUTION WITH STRESS-FREE BOUNDARIES iU(2)k k δjs =0, δjc = g ⊥ z δT . (56) ⊥ ⊥ (λ+ν(2)k2)k2 Equation (46) with eq. (42) forms a set of simulta- neous equations for T(1) and jˆ(1.5), which should be In contrast, the z component of eq. (55a) and eq. (55c) form linear homogeneous equations for δjs and δT. The supplemented by the boundary condition on jˆ(1.5). For z requirement that they have a non-trivial solution yields simplicity, we here adopt the assumption of stress-free boundaries:2 λ2+(ν(2)+κ(2))k2λ+ν(2)κ(2)(k4 R(−1)k2/k2)=0, (57) ∂2 1 − ⊥ ˆjz(1.5) = ∂z2ˆjz(1.5) =0 at z =±2. (51) withR(−1) theRayleighnumberdefinedbyeq.(35)with U′ U(2), ∆T ∆T , ν′ ν(2) and κ′ κ(2). The However, qualitative features of the convective solutions g → g → hc → → conducting solution becomes unstable when eq. (57) has will be universal among the present and more realis- a positive solution, i.e., tic/complicated boundary conditions; see the argument at the end of the preceding section. U(2)∆T (k2+k2)3 27π4 Wefirstdiscusstheheat-conductingsolutionofeq.(46) R(−1) g hc ⊥ z . (58) and its instability towards convection. We then trans- ≡ κ(2)ν(2) ≥ k2 ≥ 4 ⊥ form eq. (46) with eqs. (42) and (51) in a form suitable Thus, we have obtained the value R = 27π4/4 for to obtain periodic convective structures. c the critical Rayleigh number2,3 which corresponds to (k ,k )=(π/√2,π). ⊥ z A. Conducting solution Besidesreproducingtheestablishedresults,2theabove considerationmay also be important in the following re- spects. First, we require that: (i) terms of the z com- Let us consider the conducting solution of eq. (46) ponent of eq. (55a) all have the same order in δ with where jˆ(1.5)=0 with uniformity in the xy plane. Equa- δT =O(δ); (ii) the same be true for terms of eq. (55c). tion (46) then reduces to d2T(1)/dz2=0, which is solved This leads to the attachment of the order-of-magnitude: with eq. (42) as δj=O(δ1.5),k=O(δ−0.25)andλ=O(δ1.5),therebyjusti- 2¯j(3) fyingeq.(40d). Theconclusionk=O(δ−0.25)alsoresults T(1) =−∆Thcz, ∆Thc ≡ 3n¯κQ(2) . (52) from k = 3/2π ∼4 for the critical Rayleigh number. Second, eq. (58) removes the ambiguity in ν, κ and α to Substituting this expression into eq. (50), we obtain en- p estimate the critical Rayleigh number R . Specifically, tropy of the conducting state measured from S(0) as c we should use the mean values over 1/2 z 1/2 for S(2) = 5 (∆T )2. (53) a detailed comparison of Rc between−theory≤an≤d experi- hc −48 hc ment. 9 1 cosk C. Convective solution T1 = T˜k⊥=0,kz −k z , (64b) z Xkz We now focus on the convective solution of eq. (46) where ∆T denotes temperature difference of the con- with periodic structures. Let us introduce basic vectors hc ducting state given explicitly in eq. (52). as It follows from the energy conservation law that eq. (42b) should also hold at z = 1/2, which leads to an a =(a ,a ,0), a =(0,a ,0), a =e . (59) 1 1x 1y 2 2 3 z alternative expression for ∆T . Subtracting it from eq. cv We considerthe regioninthexy plane spannedby 1a1 (64a) yields the identity obeyed by T˜k⊥=0,kz: N and a with (j = 1,2) a huge integer, and im- N2 2 Nj T˜ k (1 cosk )=0, (65) pose the periodic boundary condition. The wave vector k⊥=0,kz z − z kisthendefinedintermsofthereciprocallatticevectors Xkz b =2π(a a )/[(a a ) a ], b =2π(a a )/[(a which is useful to check the accuracyof numericalcalcu- 1 2 3 1 2 3 2 3 1 1 × × · × × a2) a3] and b3=πez as lations. · Let us substitute eq. (61) into eq. (46). A straightfor- 3 ward calculation of using eq. (63) then leads to coupled k= ℓ b , (60) j j algebraic equations for T˜ , ˜j and ˜j as k zk pk j=1 X κ(2)k2T˜ ∆T ˜j with ℓ denoting an integer. The analysis of IVB sug- k− cv zk j § 1 gests that the stable solution satisfies |b1|∼|b2|∼π/√2. +4 T˜k′ k·j˜k′′ −δk′⊥′+k′⊥k⊥δkz′′−kz′kz Equation (56) tells us that δj⊥ and δT are out-of- k′k′′ (cid:26) (cid:20) X phase. Also noting eqs. (42) and (51), we now write δ down the steady solution of eq. (46) in the form: +δk′⊥′+k′⊥k⊥δkz′′+kz′kz + 1− k2⊥0 δk′⊥′−k′⊥k⊥δkz′′+kz′kz (cid:18) (cid:19)(cid:18) jˆ⊥(1.5)(r)= j˜⊥ksin(k⊥·r)coskzζ, (61a) −δk′⊥′−k′⊥k⊥δkz′′−kz′kz −δk′⊥−k′⊥′k⊥δkz′−kz′′kz k⊥X(6=0)Xkz (cid:19)(cid:21) + k⊥·j˜⊥k′′ −kz˜jzk′′ δk′⊥+k′⊥′k⊥δkz′−kz′′kz ˆjz(1.5)(r)= ˜jzkcos(k⊥·r)sinkzζ, (61b) (cid:16) δ (cid:17)(cid:20) k⊥X(6=0)Xkz + 1− k2⊥0 δk′⊥−k′⊥′k⊥δkz′′−kz′kz +δk′⊥′−k′⊥k⊥δkz′−kz′′kz (cid:18) (cid:19)(cid:18) T(1)(r)= T˜kcos(k⊥ r)sinkzζ ∆Tcvz T1, −δk′⊥−k′⊥′k⊥δkz′′+kz′kz =0. (66a) · − − (cid:19)(cid:21)(cid:27) Xk⊥ Xkz (61c) k2 where ζ≡z+1/2, and we have chosen ˆjz(1.5) and T(1) as ν(2)k2˜jzk−Ug(2)k⊥2T˜k even functions in the xy plane without losing the gener- 1 ality. The k summations in eq. (61) run over +4 ˜jzk′ k·j˜k′′(−δk′⊥′+k′⊥k⊥δkz′′−kz′kz k′k′′ X (cid:2) −∞≤ℓ1≤−1, 1≤ℓ2≤∞, 1≤ℓ3≤∞, (62) +δk′⊥′−k′⊥k⊥δkz′′+kz′kz +δk′⊥′+k′⊥k⊥δkz′′+kz′kz (cid:26)0≤ℓ1≤∞, 0≤ℓ2≤∞, 1≤ℓ3≤∞, −δk′⊥′−k′⊥k⊥δkz′′−kz′kz −δk′⊥−k′⊥′k⊥δkz′−kz′′kz) which covers all the independent basis functions. The +(k⊥·j˜⊥k′′ −kz˜jzk′′)(δk′⊥−k′⊥′k⊥δkz′′−kz′kz condition (39) is transformed into +δk′⊥′−k′⊥k⊥δkz′−kz′′kz +δk′⊥+k′⊥′k⊥δkz′−kz′′kz k·j˜⊥k+kz˜jzk=0. (63a) −δk′⊥−k′⊥′k⊥δkz′′+kz′kz) k Thus, j˜⊥k may be written generally as −4kz2 (k⊥·j˜⊥k′(cid:3)−kz˜jzk′)k·j˜k′′ k′k′′ j˜⊥k =−kk⊥2kz˜jzk+ ezk×k⊥˜jpk. (63b) ×(δk′⊥′+Xk′⊥k(cid:2)⊥δkz′′−kz′kz −δk′⊥′−k′⊥k⊥δkz′′+kz′kz) ⊥ ⊥ +k·j˜k′k·j˜k′′(δk′⊥′+k′⊥k⊥δkz′′+kz′kz pTehreatcuornestdainfftesre∆nTcecvbaentdweTe1ninz e=q. (16/12c)adnednottheethaveetreamge- −δk′⊥′−k′⊥k⊥δkz′′−kz′kz −δk′⊥−k′⊥′k⊥δkz′−kz′′kz) temperature shift, respectively. T±hey are fixed so as to −(k⊥·j˜⊥k′ −kz˜jzk′)(k⊥·j˜⊥k′′ −kz˜jzk′′) satisfy eq. (42) as ×(δk′⊥−k′⊥′k⊥δkz′′−kz′kz +δk′⊥′−k′⊥k⊥δkz′−kz′′kz) ∆T =∆T + T˜ k , (64a) +k·j˜k′(k⊥·j˜⊥k′′ −kz˜jzk′′)(δk′⊥+k′⊥′k⊥δkz′−kz′′kz cv hc k⊥=0,kz z Xkz −δk′⊥−k′⊥′k⊥δkz′+kz′′kz) =0, (66b) (cid:3) 10 ν(2)k2˜jpk+ 14 (ez×kk⊥)·j˜⊥k′ k·j˜k′′ Adisdafotersp:etrhioedircolslt,rtuhcetusrqeus,arweeliantvteicsetigaantde tthhee thherxeaegcoanna-l kX′k′′ ⊥ (cid:2) lattice with b1 = b2 π/√2. With these preliminaries, ×(δk′⊥′+k′⊥k⊥δkz′′−kz′kz −δk′⊥′−k′⊥k⊥δkz′′+kz′kz we trace tim|e e|vol|uti|o∼n of the expansion coefficients un- +δk′⊥′+k′⊥k⊥δkz′′+kz′kz −δk′⊥′−k′⊥k⊥δkz′′−kz′kz tiltheyallacquireconstantvalues. Choosing∆t′≤0.005 −δk′⊥−k′⊥′k⊥δkz′−kz′′kz)−(k⊥·j˜⊥k′′ −kz˜jzk′′) atinodnsℓcpr≥es5enytieedldbseeloxwce.llTenhtecinonitviaelrgsetnatceeifsorchtohseencaalscutlhae- ×(δk′⊥−k′⊥′k⊥δkz′′−kz′kz +δk′⊥′−k′⊥k⊥δkz′−kz′′kz conducting state with small fluctuations T˜k′ ∼10−2 for −δk′⊥+k′⊥′k⊥δkz′−kz′′kz +δk′⊥−k′⊥′k⊥δkz′′+kz′kz) =0, (66c) tbheeenbauspicdahtaerdmaotniecasckh.tiTmheesctoenpstbayntussi∆ngTcevqa.n(d64T).1 hAalvsoe with j˜ given by eq. (63b). Finally, entrop(cid:3)y in convec- evaluated at each time step is entropy measured with ⊥k respect to the heat-conducting state: tion is obtained from eqs. (50) and (61c) as ∆S S(2) S(2), (68) 5 (∆T )2 1 ≡ cv − hc S(2) = cv T2+ T˜ 2(1+δ ) cv −4" 12 − 1 4 k | k| k⊥0 where Sh(2c) and Sc(v2) are given by eqs. (53) and (67), re- X spectively. Wetherebytracetimeevolutionof∆S simul- 1+cosk +∆T T˜ z . (67) taneously. The above procedure is carried out for each cv k⊥=0,kz kz # fixed periodic structure. Xkz We have studied the range: 1 R(−1)/R 10. c ≤ ≤ Although the region extends well beyond the Busse V. NUMERICAL RESULTS balloon3,6 of stability for classical gases,5,8,9 it will be worth clarifying the basic features of steady periodic solutions over a wide range of the Rayleigh number. We now present numerical results on entropy of Rayleigh-Bénard convection obtained by solving eqs. (66a)-(66c). As a model system we consider Ar at B. Results T¯=273Kunderatmosphericpressureandusethe values (37) for the parameters in eqs. (66a)-(66c). We also fix the heat flux density ¯j at z = d/2 so that tempera- Figure 1 plots ∆S as a function of t′ for the Rayleigh Q − number R(−1)=1.2R which is slightly above the criti- ture difference ∆T of the conducting state, expressed c in terms of ¯j as hecq. (52), takes the value ∆T =1K. cal value Rc=27π4/4. The letters r, s and h denote (r) Q hc The Rayleighnumber R(−1) is then controlledby chang- ing the thickness d. We here adopt the units described -7 x10 above eq. (32) with k =1. B Recent experiments on Rayleigh-Bénard convection 3 havebeenperformedmostactivelywithcompressedclas- sical gases5,8,9 where controlled optical measurements of r s convectivepatternsarepossible. Ithas beenpointedout 2 h thatthebasicn( l−1)andT dependencesofeq.(29)are well satisfied even∝for those gases under high pressures.8 S Thus, the present consideration with eq. (37) has direct D relevance to those experiments. 1 (-1) R =1.2R c A. Numerical procedures 0 0 20 40 60 80 100 To solve the coupled equations numerically, we first t' multiply eqs. (66a)-(66c) by 108, 1010 and 1010, respectively, and rewrite them in terms of T˜′ 103T˜ and k ≡ k ˜j′ 105˜j to obtain equations of O(1). We next replace zke∂r≡o˜js′ o/f∂ttkh′e,rreigsphte-chtaivnedlys.idTehsebsye−tim∂T˜ek′d/e∂rti′v,a−ti∂v˜jezs′ka/r∂etw′ahnadt FtoIGth.e1:heaTti-mcoendevuocltuintigonstaotfeefnotrroRp(y−1m)e=as1u.2reRdc.wTithherleestpteecrst − pk r, s and h denote roll, square and hexagonal, respectively, onewouldgetbyretainingtimedependenceintheexpan- distinguishinginitialfluctuationsaroundtheheat-conducting sion coefficients of eq. (61). We then discretize the time solution; see text for details. The final state of t′ & 80 is derivatives as ∂T˜′(t′)/∂t′ [T˜′(t′+∆t′) T˜′(t′)]/∆t′, the roll convection, whereas the intermediate plateaus of s for example. Theksummat≈ionskover k are−trukncated by and h correspond to the square and hexagonal convections, using a finite value ℓ (&5) in place of in eq. (62). respectively. c ∞