Table Of Content

Exactlinearhydrodynamics from the Boltzmannequation I.V. Karlin,1,∗ M. Colangeli,2 and M. Kro¨ger2,† 1AerothermochemistryandCombustionSystemsLab, ETHZu¨rich, CH-8092Zu¨rich, Switzerlandand School ofEngineeringSciences, UniversityofSouthampton, SO171BJ,Southampton, UnitedKingdom 2PolymerPhysics,DepartmentofMaterialsandMaterialsResearchCenter,ETHZu¨rich,CH-8093Zu¨rich,Switzerland (Dated: 2008-01-1809:47:09) Exact(toallordersinKnudsennumber)equationsoflinearhydrodynamicsarederivedfromtheBoltzmann 8 kineticequationwiththeBhatnagar-Gross-Krookcollisionintegral.Theexacthydrodynamicequationsarecast 0 inaformwhichallowsustoimmediatelyprovetheirhyperbolicity,stability,andexistenceofanH-theorem. 0 2 PACSnumbers: 51.10.+y(Kinetictheory)05.20.Dd(Kinetictheory) n Keywords:Kinetictheory;heattransfer;hydrodynamics;hyperbolicequations;Boltzmannequation a J Hydrodynamicsassumesthatastateofafluidissolelyde- 8 0 1 scribed by five fields: density, momentum and temperature. DerivationoftheNavier-Stokes-Fourier(NSF)hydrodynamic −0.2 Re(ω ) ] equations from the Boltzmann kinetic equation as the first- ac h Im(ω ) / 10 c order approximation in the Knudsen number (ratio between −0.4 Re(ωac) hydrodynamic modes e mean free path and a flow scale) by Enskog and Chapman −0.6 Im(ωaacc) / 10 (real and imaginary parts) m is a textbook example of a success of statistical physics [1]. Re(ω ) diff - Recentrenewedinteresttotheproblemsbeyondthestandard −0.8 Im(ωdiff) t hydrodynamicsisdue,inparticular,toflowsimulationandex- Re(ω ) a shear t perimentsatamicro-andnano-scale[2–4]. However,almost −1 Im(ωshear) .s acenturyofefforttoextendthehydrodynamicdescriptionbe- 0 0.5 1 k 1.5 2 2.5 t a yondtheNSFapproximationfailedeveninthecaseofsmall m deviationsaroundthe equilibrium. Inordertoappreciatethe FIG. 1: Exact hydrodynamic modes ω of the Boltzmann-BGK kinetic problem, let us remind that, in the NSF approximations, the - equation asafunction ofwavenumberk (twocomplex-conjugated acous- d decay rate of the hydrodynamic modes is quadratic in the ticmodesωac,twicedegeneratedshearmodeωshearandthermaldiffusion on wavevector,Re(ω) ∼ −k2,andisunbounded. Ontheother modeωdiff).Thenon-positivedecayratesRe(ω)attainthelimitofcollision hand, Boltzmann’s collision term features equilibration with frequency(−1)ask→∞. c finite characteristicrates. This“finite collisionfrequency”is [ obviously incompatible with the arbitrary decay rates in the 1 NSFapproximation:intuitively,hydrodynamicmodesatlarge v popular in applications [9], and features a single relaxation k cannot relax faster than the collision frequency. Now, the 2 rate. Theresultforthehydrodynamicmodesisdemonstrated classicalmethodofEnskogandChapmanextendsthehydro- 3 inFig.1. ItisclearfromFig.1thattherelaxationofnoneof dynamicsbeyondtheNSFinsuchawaythatthedecayrateof 9 thehydrodynamicmodesisfasterthanω = 1whichisthe 2 thenextorderapproximations(Burnettandsuper-Burnett)are − collision frequencyin the units adopted in this paper. Thus, . polynomialsofhigherorderink. Insuchanextension,relax- 1 theresultfortheexacthydrodynamicsindeedcorrespondsto ation rate may becomecompletelyunphysical(amplification 0 theaboveintuitivepicture.Below,weapplythemethodofin- 8 instead of attenuation), as first shown by Bobylev [5] for a variantmanifold [10] to derive the hydrodynamicequations. 0 particularcaseofMaxwellmolecules.Thisindicatesinability Thenon-perturbativederivationismadepossiblewithanop- : oftheChapman-Enskogmethodtotackletheaboveproblem, v timalcombinationof analyticalandnumericalapproachesto andnon-perturbativeapproachesaresought. Theproblemof i X exact hydrodynamics has been studied in depth recently for solvetheinvarianceequation. PointofdepartureisthelinearizedBoltzmann-BGKequa- r toy(finite-dimensional)models-momentsystemsofGrad-in a tionforthedeviation f =f fGMofthedistributionfunc- [6–8],andmanyremarkableresultswereobtained. Inpartic- tionf fromaglobalM△axwelli−anfGM(c2) = π−3/2e−c2. In ular,in[6,7]itwasshownthattheexacthydrodynamicequa- thereciprocalspace,itreads, tions are hyperbolicand stable for all wave numbers. How- ever,for“true”kineticequationssuchquestionsremainopen. ∂ f = ik c f δf; δf =f fLM, (1) t △ − · △ − − InthisLetterwederiveexacthydrodynamicequationsfrom withthewavevectork = ke defininge ,k k,peculiar thelinearizedBoltzmannequationwiththeBhatnagar-Gross- k k ≡ | | velocity c and time t. All quantities are considered dimen- Krook (BGK) collision term. This kinetic equation remains sionless, i.e.,reducedwiththeunitsoftherelaxationtimeτ, the thermalvelocity 2kBT/mand mass m of the particle. In(1),thelinearizedlpocalMaxwellianisfLM =fGM(1+ϕ0) ∗†UElRecLt:rownwicw.adcdormespsl:[email protected] ewrhagerees1ar+edϕe0fi=nehd1fiofr+ar2bcitr·ahrcyifω+vi32a(cω2−f 23)hc2ω−f32di3fc.,Aanvd- h i ≡ R 2 weintroducepertinentquantitieswhichcharacterizedeviation σ1k σ2k σ3k σ4 fromtheglobalequilibrium:n 1 1(densityperturba- tion),u ≡ hcif (velocitypertur≡bahtioinf)−andT ≡ 32hc2− 23if hλ−kkδ2XB1i hλikkδAX2i hλ−kkδ2XC3i hckick⊥DδY4i (temperatureperturbation). Sincethescalarproductbetween kandcappearsin(1),thedistributionfunctionofferssymme- B0 =−31 A0 =−23 C0 =0 D0 =−21 trywithrespecttothee -axis,whichisnotuniaxialincaseu real,⊕ imag,⊕ real,⊕ imag,⊖ k is not collinear with ek. We denote the two components of q1k q2k q3k q4 the mean velocity as uk = u·ek and u⊥ = u·e⊥, where hγkδX1i hγkδX2i hγkckδX3i h(c2− 25)c⊥δY4i theunitvectore⊥belongstotheintersectionoftheplaneper- ikX −k2Z ikY −k2U pue=ndiuckuelar+tou⊥kean.dEqthueatpiolannseofspcahnannegdefboyrmkoamndenuts, sωo tharaet X0 =0 Z0 =−16 Y0 =−45 U0 = 21 obtainedkbyinteg⊥rationoftheweighted(1)overd3cash i imag,⊖ real,⊖ imag,⊖ real,⊕ ∂ ω = ik cω ω . (2) TABLEI:Symmetryadaptedcomponentsof(nonequilibrium)stress th i△f − ·h i△f −h iδf tensor σ andheat fluxq, introducedin(7a)and(7b), respectively. In order to calculate such averages, we can switch to spher- Row2:Microscopicexpressionofthesecomponents(averagingwith ical coordinates. For each (at present arbitrary) wave vec- theglobalMaxwellian).Short-handnotationused:λk =c2−c2 and k 3 tor, we choose the coordinate system in such a way that its γk = (c2− 5)c . Row3: Expressionofthecomponents interms 2 k z-direction aligns with e . We can then express c in terms of(asweshow,real-valued)functionsA–Z (seetext). Row4: Val- k of its norm c, a vertical variable z and plane vector e (az- uesoffunctionsA–Zatk=0.Thesevaluesrecoverhydrodynamic φ imuthal angle e e = cosφ) for the present purpose as equationsuptoBurnettapproximation. Row5: Paritywithrespect c/c = √1 z2φe· +⊥ze . We could have equally chosena toz–symmetric(⊕)orantisymmetric(⊖)–ofthepartofthecor- − φ k responding δX entering the averaging in row 2, and whether this fixedcoordinatesystemintheplaneorthogonaltok,andtwo partisimaginaryorreal-valued(seeFig.2). Row3isanimmediate fieldsinsteadofu⊥plusanangle,viz.u⊥+iu⊥ =eiφu⊥.Due x y consequenceofrow5. to isotropy, u⊥ alone fully represents the twice degenerated (shear)dynamics. Inordertosimplifynotationandcompute thedynamicsofallfivefieldsweintroduceafour-dimensional whichisanonlinearintegralequationfortheunknownfields v(sxiemckt,poxlre4)foowfrmihth,yϕdxr0ko=d≡yXn(an0m,uixck.,fiTTeh)ledasvn,edcxtxo4r≡X=0u(c⊥xa.n1T,ixmh2em,nxeϕ3d0,iaxtta4ek)lyesb=ae δoHXfe(r,2e,)b,eωfcoairussaseismu∂ititlxaabrhlwyasicthhtooXsbe0enarvenepdcltadocireffωderbsfyuflrfitohlmleinrXgigh0hωtmih△aainfndly=sbidxee-. readoff,wehaveX0(c·,z)=(1,2c ,c2 3,2c ),wherewe causeofconventionsforprefactorsinthetemperatureandve- introduced,forlateruse,theabbrevikation−s 2 φ locitydefinitions,ω =(1,c ,2(c2 3),c ). Withinthesame k 3 −2 φ c eigen-closure,Eq.(2)islinearinxandhencewrittenas c c e , c c e , c φ , (3) k ≡ · k φ ≡ · ⊥ ⊥ ≡ e⊥·eφ ∂tx=M x. (6) · suchthatik c = ikc , c = c e +c e with c = cz and · k ⊥ φ k k k The matrix M solely depends on the non-hydrodynamic c⊥ = c√1−z2, contrasted by cφ (and eφ), do not depend fields, the heat flux q c(c2 52) f and the stress tensor ontheazimuthalangle. Similarly,weintroduceyetunknown ≡ h − i σ cc , where s denotesthe symmetric tracelesspart tfiheelddsisδtXrib(uct,ikon)wfuhniccthiocnh,aδraϕct=erizδeft/hfeGnMonienquteirlimbrsiuomf tphaerthoyf- of≡a tehnsorisf[6, 11], s = 12(s+sT)− 13tr(s)I. Using(4), thestresstensorandheatfluxuniquelydecomposeasfollows drodynamicfieldsxthemselves, δϕ=δX x=δX1n+δX2uk+δX3T +δX4u⊥, (4) 3 · σ = σk e e +σ⊥2e e , (7a) Twhhiesre“eδiXge4n”f-accltoosruizrees(a4s) δwXhi4c(hc,fzo,rφm)al=ly a2nδdY4v(ecr,yz)geenφe·rael⊥ly. q = qke2 +kqk⊥e , k ⊥ (7b) k ⊥ addressesthefact,thatwewishtonotincludeotherthanhy- drodynamicvariablesimpliesa closurebetweenmomentsof withthemomentsσk =(σ1k,σ2k,σ3k)·xkandσ⊥ =σ4x4,and similarlyforq(seeRow2ofTab.I).Theprefactorsarisefrom thedistributionfunction,tobeworkedoutindetailbelow. It theidentities e e : e e = 2 and e e : e e = 1. The assumes the existence of an invariant manifold, and the hy- k k k k 3 k k k ⊥ 2 appearanceof δY4 rather than δX4 in the expression for the drodynamic fields as slow variables which leave the higher orthogonalmoment(inTab.I)reflectsthefactthatwehaveal- moments “slaved”. In order to ensure these contributionsto readyintegratedouttheangularvariable, 2πe e e dφ= notinterferewiththelocalMaxwellian,onehasthefreedom R0 φ φ· ⊥ torequire X0 δf = 0withoutproducinganylimitation,i.e., πe⊥. We note in passing that, while the stress tensorhas, in h i general, threedifferenteigenvalues,in the presentsymmetry bykeepingxtobedefinedthroughthelocalMaxwellianpart adaptedcoordinatesystemitexhibitsavanishingfirstnormal of the distribution function. Using the above form for δf in (1), and using the canonicalabbreviation X X0 +δX, stressdifference. Sincetheintegralkernelsofallmomentsin △ ≡ (7) do not depend on the azimuthalangle, these are actually yields two-dimensionalintegrals over c [0, ] and z [ 1,1], X ∂ x= ikc X x δX x, (5) weightedby2πc2fGM(c2)δX . ∈ ∞ ∈ − △ · t − k△ · − · µ 3 Stresstensorandheatfluxcanyetbewritteninanalterna- [6,7]andbelow),oncethefunctionsA–Z areexplicitlyeval- tive form which is defined by Row 3 of Tab. I. As we will uated. For that, we still requireδX. Combining(5)and (6), seelateron,duetobasicsymmetryconsiderations,thehereby andrequiringthatthe resultholdsforanyx(invariancecon- introducedfunctions A–Z are real-valued. We postpone the dition),weobtainaclosed,singularintegralequation(invari- related proof, and proceed by using these functions A–Z to anceequation)forcomplex-valuedδX, splitMintopartsasM=Re(M) iIm(M)with − δX=X0 (M+[ikc +1]I)−1 X0. (9) · k − 0 0 0 0   Notice that δX vanishes for k = 0, and that (9) is supple- Re(M) = k2 320X A0 230Y 00 , (8) mvaennistehdinwgitLhagthraenbgaesimcuclotinpsltirearisnt(mhXat0riixδf) =X00δ,Xor,ewquhailclhy,,   h i  0 0 0 D  however, is automatically dealt with if we only evaluate anisotropic (irreducible) moments with δf, such as those 0 1 0 0   listed in Tab. I. The implicit equation (9) is identical with Im(M) = k 12−0k2B 32(1−0k2Z) 12−0k2C 00 . trheseuelti,gwenit-hclMosufrreo(m4)(,8a)n,dcki,sXo0u,ramndaiAn–aZnddperfiancetidcainllyanudsejufustl  0 0 0 0  before(3)andTab.I,respectively. Weiterativelycalculate(i)δXdirectlyfrom(9)foreachk intermsofM,(ii)subsequentlycalculatemomentsfromδX bynumericalintegration. Importantly,thefix pointoftheit- eration(i)-(ii)-(i)-..isuniqueforeachk,i.e.,doesnotdepend ontheinitialvaluesformomentsA–Z,aslongaswechoose real-valued ones which are consistent with (9), as we prove in the next paragraph. In addition, two other computational strategies were implemented: First, we used continuation of functionsA–Z fromtheirvaluesatk =0tosolve(9)withan incrementalincreaseofk,wherethesolutionatkwasusedas theinitialguessfork+dk. Second,weusedalsoacontinu- ation“backwards”inwhichthesolutionatsomek (obtained byconvergentiterationswitharandominitialcondition)was usedastheinitialguessfora solutionatk dk. Both these − strategiesreturnedthesamevaluesoffunctionsA–Z ascom- putedbyiterationsfromarbitraryinitialcondition. Thesolu- tionδXallowstocalculatethewholedistributionfunctionf via (4) as illustrated by Fig. 2. For resulting moments for a widerangeofk-valuesseeFig.3. Finally, we need to clarify the origin of row 5 in Tab. I (whichisdirectlyillustratedbyFig.2)anditsimplicationon thestructureofM(8)whoseentriesare–apriori–complex- valued functions to be calculated with complex-valued δX. Wewishtomakeuseofthefactthatallintegralsoverz van- ish for odd integrands. To this end we introduce abbrevia- tions ( ) for a real-valued quantity which is even (odd) ⊕ ⊖ with respect to the transformation z z. One notices X0 = ( , , , ), and we recall th→at A−–Z are integrals ⊕ ⊖ ⊕ ⊕ overeitherevenoroddfunctionsinz,timesa componentof δX (see Tab. I). Let us prove the consistency of the specified symmetry of M and the invariance condition: Start by assumingA–Z to be real-valuedfunctions. ThenM = µν ⊕ if µ + ν is even, and M = i otherwise. This implies FIG.2: (Coloronline)Sampledistributionfunctionf(c,k)atk=1,fully µν ⊕ characterizedbythefourquantitiesδX1,2,3(c,z)andδY4(c,z).Shownhere δX1 = ⊕ + i⊖, δX2 = ⊖ + i⊕, δX3 = ⊕ + i⊖, and areboththeirreal(left)andimaginaryparts(rightcolumn). Inordertoim- δX4 = +i , i.e., different symmetry properties for real ⊕ ⊖ provecontrast,weactuallyplotln|1+fGMδXµ|multipliedbythesignof andimaginaryparts. Withthese“symmetry”expressionsfor δXµ.Samecolorcodeforallplots,rangingfrom−0.2(red)to+0.2(blue). X0,δX,andMathand,wecaninsertintotherighthandside of the equation, δX = (X0 +δX) (M+i I), which is · ⊖ identicalwiththeinvarianceequation(9). Thereareonlytwo Note that the checkerboard structure of the matrix M (8) cases to consider, because M has a checkerboard structure, is particularlyusefulfor studyingpropertiesof the hydrody- i.e.,onlytwotypesofcolumns: Columnsµ = 1andµ = 3: namic equations(6), such as hyperbolicityand stability (see δXµ = +i becauseM1−3,4 = 0; Columnsµ 2,4 : ⊕ ⊖ ∈ { } 4 δXµ = +i ifMµ,1−3 =0(whichisthecaseforcolumn ∂tx′ =M′ x′,withM′ =T M T−1ismanifestlyhyper- 4)and ⊕+i ⊖ifMµ,4 =0(whichisthecaseforcolumn2). bolicandsta·ble;Im(M′)issym· me·tric,Re(M′)issymmetric ⊖ ⊕ and non-positive semi-definite. The corresponding transfor- mationmatrixTcanbeeasilyreadofftheresultsobtainedin 0.5 [6,7]forGrad’ssystemssincethestructureofthematrixM (8) is identical to the one studied in [6, 7]. We have explic- itlyverifiedthatmatrixT(equations(21)–(23)inRef.[7]and 0 Z o (13)inRef.[6])withthefunctionsA–Zderivedhereinisreal- A t valued and thus render the transformedhydrodynamicequa- nts −0.5 A tionsmanifestlyhyperbolicandstable.Wenotethatthisresult e B m – hyperbolicity of exact hydrodynamic equations – strongly C o m D supportsarecentsuggestionbyBobylevtoconsiderahyper- −1 U X bolicregularizationof the Burnettapproximation[12]. Sim- Y ilarly, using the hyperbolicity, an H-theorem is elementary Z −1.5 provenasin[6,12]. Finally,usingtheaccuratedataforfunc- 0 0.5 1 1.5 2 2.5 3 k tionsA–Z, wecanwriteanalyticapproximationsforthehy- drodynamicequations(6)insuchawaythathyperbolicityand FIG.3:MomentsA–Zvs.wavenumberkobtainedwiththesolutionof(9). stabilityisnotdestroyedinsuchanapproximation(see[7]). In conclusion, we derived exact hydrodynamic equations Wehavethusshownthatbothsidesoftheinvarianceequa- fromthelinearizedBoltzmann-BGKequation.Themainnov- tion(9)haveequalsymmetryproperties,andthatδXwiththe elty is the numericalnon-perturbativeprocedureto solve the specified symmetriesis consistentwith real-valuedmoments invariance equation. In turn, the highly efficient numerical A–Z. The proof implies, that any iteratively obtained solu- approach is made possible by choosing a convenient coor- tion, if it exists, starting with arbitrary real-valued moments dinate system and establishing symmetries of the invariance A–Z in (9) to evaluate δX must convergeto real-valuedso- equation. Theinvariantmanifoldin the spaceof distribution lution A–Z. Since the solution is smoothly varying with k, functionsistherebycompletelycharacterized,thatis,notonly and since A–Z at k = 0 are known and are real-valued, the equationsof hydrodynamicsare obtained but also the corre- momentsmustbereal-valuedoverthewholek-space. sponding distribution function is made available. The perti- WiththeresultforthefunctionsA–Z athand,theextended nentdatacanbeused,inparticular,asamuchneededbench- hydrodynamic equations are closed. Let us briefly discuss markforcomputation-orientedkinetictheoriessuchaslattice the pertinent properties of this system. First, the general- Boltzmannmodels, as well as forconstructingnovelmodels ized transport coefficients are given by the nontrivial eigen- usingquadraturesin the velocityspace [13, 14]. Finally, we values of k−2Re(M): λ2 = A (elongation viscosity), have established a novel non-perturbativecomputationalap- λ3 = −32Y−(thermal diffusivity),−and λ4 = −D (shear vis- proachtofindinginvariantmanifoldsofkineticequations. cosity). All these generalized transportcoefficients are non- negative (see Fig. 3). Second, computing the eigen-values The present approach can be extended to the Boltzmann of matrix M we obtain the dispersion relation ω(k) of the equationwith othercollision terms. Theabovederivationof correspondinghydrodynamicmodesalreadypresentedinFig. hydrodynamicsisdoneunderthe standardassumptionof lo- 1. Third, a suitable transform of the hydrodynamic fields, calequilibrium[1], howeverthe assumptionitself is open to x′ = T x, where T is a real-valued matrix, can be es- furtherstudy[15][WethankH.C.O¨ttingerforthisimportant · tablishedsuchthatthetransformedhydrodynamicequations, remark].I.V.K.acknowledgessupportofCCEM-CH. [1] S. Chapman, T. G. Cowling, The Mathematical Theory of [8] A.N.Gorban,I.V.Karlin,Phys.Rev.Lett.77(1996)282. NonuniformGases (CambridgeUniv.Press,NewYork,1970). [9] C.Cercignani,TheoryandApplicationoftheBoltzmannEqua- [2] A.Beskok, G.E.Karniadakis, Microflows: Fundamentalsand tion (ScottishAcademicPress,Edinburgh,1975). Simulation (Springer,Berlin,2001). [10] A.N.Gorban,I.V.Karlin,Transp.Th.Stat.Phys.23(1994)559. [3] S. Ansumali, I.V. Karlin, S. Arcidiacono, A. Abbas, N.I. [11] M. Kro¨ger, Models for Polymeric and Anisotropic Liquids Prasianakis,Phys.Rev.Lett.98(2007)124502. (Springer,Berlin,2005). [4] D.M.Karabacak,V.Yakhot,K.L.Ekinci,Phys.Rev.Lett.98 [12] A.V.Bobylev,J.Stat.Phys.124(2006)371. (2007)254505. [13] S.S. Chikatamarla, I.V. Karlin, Phys. Rev. Lett. 97 (2006) [5] A.V.Bobylev,Sov.Phys.Dokl.27(1982)29. 190601. [6] M.Colangeli, I.V.Karlin, M.Kro¨ger, Phys.Rev.E76(2007) [14] X.Shan,X.He,Phys.Rev.Lett.80(1998)65. 022201. [15] H.C. O¨ttinger, H. Struchtrup, Multiscale Model. Simul. 6 [7] M.Colangeli, I.V.Karlin, M.Kro¨ger, Phys.Rev.E75(2007) (2007)53. 051204.