Printed3January2013 Molecular dissipation in the nonlinear eddy viscosity in the Navier-Stokes equations: modelling of accretion discs 3 G. Lanzafame⋆ 1 0 INAF - Osservatorio Astrofisico di Catania, ViaS. Sofia 78 - 95123 Catania, Italy 2 n a J Accepted ——-.Received——-;inoriginalform——- 1 ] ABSTRACT n Physicaldamping,regardingthenonlinearNavier-Stokesviscousflowdynamics,refers y d to a tensorial turbulent dissipation term, attributed to adjacent moving macroscopic - flow components. Mutual dissipation among these parts of fluid is described by a u braking term in the momentum equation together with a heating term in the energy l f equation, both responsible of the damping of the momentum variation and of the s. viscous conversion of mechanical energy into heat. c A macroscopic mixing scale length is currently the only characteristic length i needed in the nonlinear modelling of viscous fluid dynamics describing the nonlin- s y ear eddy viscosity through the kinematic viscosity coefficient in the viscous stress h tensor, without any reference to the chemical composition and to the atomic dimen- p sions.Therefore,inthis paper,we write a new formulationforthe kinematic viscosity [ coefficienttotheturbulentviscousphysicaldissipationintheNavier-Stokesequations, 3 where molecular parameters are also included. v Resultsof2Dtestsareshown,wherecomparisonsamongflowstructuresaremade 8 on 2D shockless radial viscous transport and on 2D damping of collisional chaotic 4 turbulence.Anapplicationtothe3Daccretiondiscmodellinginlowmasscataclysmic 8 variables is also discussed. 6 Consequencesofthekinematicviscositycoefficientreformulationinamorestrictly 0. physical terms on the thermal conductivity coefficient for dilute gases are also dis- 1 cussed. 2 The physical nature of the discussion here reported excludes any dependence by 1 the pure mathematical aspect of the numerical modelling. : v Key words: accretion, accretion discs – hydrodynamics – binaries: close – stars: i X novae, cataclysmic variables. r a 1 INTRODUCTION their contemporary heating. A turbulent kinematic viscos- itycoefficientischaracterizedbyamacroscopicscalelength Physicaldissipationintheviscousfluiddynamicsistheonly multiplied by a scale velocity in the Von K´arma´n descrip- physical mechanism for the attenuation of the momentum tion, originally formulated to describe a repeating pattern transferandfortheconversionofmechanicalenergy(kinetic of swirling vortices caused by the unsteady separation of +potential)intoheat.Itoriginates from microscopic parti- flow of a fluid over bluff bodies. Hence, a mixing length is cle interactions on molecular scale lengths. For this reason, often required in the formulation of the kinematic viscos- it is not currently included in the nonlinear Navier-Stokes ity coefficient ν in the viscous stress tensor describing the equationsfor viscous flows, where macroscopic spatial reso- nonlinearturbulenteddyviscosity,withoutanyreferenceto lution lengths of moving fluid components are much larger thechemical composition and to the microscopic molecular than molecular scale lengths. dimensions.Moreover,physicalturbulentviscosity often in- Thus, a physical turbulent viscosity is used in the cludes arbitrary parameters, to be set case by case, as it is Navier-Stokes equations, as a tensorial viscous dissipation for of the well known Shakura (1972); Shakura& Sunyaev term, relative to mutual interactions among contiguous (1973) formulation for disc structures. macroscopicmovingflowparts,producingtheirbrakingand Adissipationmechanismisalwaysnecessaryinthecom- putationalcollisionalfluiddynamics,eveninthenonviscous ⋆ E-mail:[email protected] modellingtosolvethestrictlyhyperbolicEulerequations,if 2 G. Lanzafame flowdiscontinuities(theRiemannproblem) mustbesolved. 2 DISSIPATION IN VISCOUS AND NON In the physically inviscid fluid dynamics, ”shock captur- VISCOUS FLUID DYNAMICS ing” methods adopt either an artificial viscosity contribu- In the physically non viscous flows, the hyperbolic Euler tion or take advantage of some moderate intrinsic numeri- system of equations cal dissipation. Instead, ”shock tracking” methods develop a dissipation separately handling shock fronts using appro- dρ +ρ∇·v=0 continuity equation(1) priate Riemann solver algorithms, through algebraic aver- dt ages between the two left-right sides as Rankine-Hugoniot dv ∇p =− +f momentum equation(2) jump conditions, in the Godunov-type methods. In addi- dt ρ tion, a further dissipation is unavoidably intrinsically gen- d 1 1 erated by truncation errors (Fletcher 1991; Hirsch 1997). ǫ+ v2 =− ∇·(pv)+f ·v energy equation(3) dt 2 ρ In the finite difference methods, dissipation comes about (cid:16) (cid:17) by second order derivatives coming from the Taylor series dr =v kinematic equation.(4) expansion for incremental ratios of the first order spatial dt derivatives (Park & Kwon 2003), especially for implicit in- mustbesolved,togetherwith thestateequation(EoS) tegration techniques.Suchadissipation contributionisalso of thefluid usefultosmoothoutspuriousheatingandtotreattransport p=f(γ,ρ,ǫ,r,v) state equation(5) phenomena. Most of the adopted symbols have the usual meaning: Thephysical,turbulentdissipationandthefictitiousar- d/dt stands for the Lagrangian derivative,ρ is thegas den- tificialornumericalonesareconceptuallydistinct,although sity, ǫ is the thermal energy per unit mass, p is the ideal formally similar, as shown in Molteni et al. (1991); Murray gas pressure, here generally expressed as a function of local (1996); Okazaki et al. (2002), and discussed by Lanzafame properties,vandrarethevectorvelocityandposition,f is (2008, 2009) in the case of smooth particle hydrodynamics theexternalforcefieldperunitmass.Theadiabaticindexγ (SPH)accretiondiscmodellinginclosebinaries(CBs),orby has the meaning of a numerical parameter whose value lies Toro (1999) in the case of finite difference or for finite vol- in the range between 1 and 5/3. ume integration techniques. In both cases, some arbitrary SincetheRiemannproblemmustbecorrectlysolvedfor parameters,tobetunedcasebycase,areincludedand/ora collisional flows in the case of shocks, a dissipation mecha- dependence on the spatial resolution length affects the nu- nism is necessary otherwise frontal colliding flows trespass merical dissipation. each other. Such a dissipation could be either explicit, as an artificial viscosity term for shock capturing schemes, es- pecially for finite volume schemes, or it could be intrin- Therefore, in this paper we propose the formulation of sic through a specific Riemann solver code (LeVeque 1992; a physical turbulent kinematic viscosity coefficient in the Fletcher 1991; Hirsch 1997; LeVeque 2002; Park & Kwon Navier-Stokes equations, free of any arbitrary parameters, 2003)eitherforshocktrackingschemesorforEulerianfinite wheremicroscopicphysicalcharacteristicsarealsoincluded. difference schemes by commuting mathematical derivatives Amacroscopicscalelength(mixinglength)isobviouslystill inincrementalratios.Asanexample,inthefinitedifference used, since solutions of the nonlinear Navier-Stokes equa- techniques,theconversionofthe1storderspatialderivative tions, involve macroscopic physical properties. At the same time, a reformulation for ν in a more physical sense is also of the generic physical quantityu is (ui−ui−1)/∆x, where iisaspatialgridindexand∆xisthegridspatialresolution coherenttoareformulationforthethermalconductivityco- length.Forabetterstableresult,thesameincrementalratio efficient c for dilute gases. is rewritten as In §2 of this paper we discuss some general aspects of ui−ui−1 = ui+1−ui−1 − ui+1−2ui+ui−1, (6) ∆x 2∆x 2∆x dissipation in the computational collisional fluid dynamics; which gives a much better stability, at the cost of a in §3, we shortly describe characteristics of some adopted reduction in accuracy. The contribution of the second term physical turbulent kinematic viscosity coefficients ν, while in§4weformulatethephysicaldevelopmentofbothν andc (ui+1−2ui+ui−1)/∆xanalyticallycorrespondstoa2ndor- derspatialderivative,workingexactlylikearealviscousnon coefficientsincludingmicroscopicmolecularcharacters.In§5 physicalcontribution.Thesemanipulationsofspatialderiva- weshowsomeresultsforsomeessential2Dtestsonshockless tives in the incremental ratios in finite terms are necessary radialviscoustransportinanannularring,aswellasforthe to ensure stability to the solutions of hyperbolic systems of dampingof2DBurger’sturbulence.Instead,in§6,weshow equations. This means that the inclusion, or the numerical an astrophysical application in the case of a 3D accretion developmentofanonphysicaldissipationdistortsnumerical disc modelling in a low mass close binary (LMCB) system, results especially for non collisional events like shear flows comparing 3D accretion disc structures obtained by using or transport phenomena. These numerical difficulties arise different ν coefficients and gas compressibility. when the EoS Despite the adoption of a numerical techniqueintrinsi- p=(γ−1)ρǫ perfect gas equation(7) cally quite viscous, like the SPH, successful physically vis- is adopted for ideal flows. Instead, an EoS as: cous results unconditionally show that thenew formulation forthekinematicviscositycoefficientworkswellwithoutany p∗ = ρc2 1−Cn−1/3∇·v 2, (8) restriction. γ s 3c s (cid:18) (cid:19) Molecular dissipation and macroscopic physical viscosity 3 includes a real macroscopic physical dissipation cor- 3 THE PHYSICAL KINEMATIC VISCOSITY rectly handling the Riemann problem, as well as transport COEFFICIENT andshearflowsfreeofanylocalgascompression(Lanzafame Typical kinematic laboratory viscosities are of the order of 2010a,b; Lanzafame et al. 2011). c is the sound velocity, n s ν =0.001−1 cm2 s−1, to be compared with inertial forces is thenumerical density,while in the ratio C = 1arccot DvR , (9) Re= inertial forces ≡ lflowvflow (14) π cs viscous forces ν (cid:16) (cid:17) whereD≫1andwherev isthecomponentofvelocity where l and v arethecharacteristic length and R flow flow alongthedirectionofcollision.Disalargenumberdescrib- velocity scales of the microscopic flow. Laboratory experi- inghowmuchtheflowdescriptioncorrespondstothatofan ence shows that for Re>Re ≈102−103, flow becomes crit ideal gas: D ≈ λ/d, being λ ∝ ρ−1/3 the molecular mean turbulent. Re is the characteristic Reynolds number as crit free path, and being d the mean linear dimension of gas observed so far. molecules. The physical dissipation, expressed by the two In the full nonlinear approach, the full non linearity furtherterms in eq.(8) (thelinear and thequadraticterms of the Navier-Stokes equations is considered, where spatial in∇·v)ofthereformulated EoS,bettertreatsbothshocks derivatives of the entire velocity field are used. Neither the and shear flows, even in a Lagrangian description. Their Reynolds averages of the Navier-Stokes equations in boxes inclusion substitutes artificial viscosity terms and does not of intermediatesize (as in thelinear approach), northefull representa physicalturbulentviscouscontribution,butthe Navier-Stokes equations working with spatial gradients of real physical dissipation coming out because eq. (7) of the the mean velocity field (as in the nonlinear Boussinesq ap- EoS should strictly be applied only to macroscopic static proach (Schmitt 2007)), are considered. or quasi-static processes. Notice that in this case, this real To characterize a nonlinear macroscopic physical kine- macroscopic physical dissipation does not originate from a maticviscositycoefficientν,characteristiclengthandveloc- physicalviscosity.Instead,itoriginatesfromtheirreversible ityscales l andv areneeded,whichareunknown,inprinci- thermodynamic process and is better evidenced in a La- ple. grangian description. Typically(Prandtl1925),amixinglengthmodelcanbe Inthephysicallyviscousflows,theNavier-Stokesequa- used,where tions explicitly include macroscopic physical dissipation ∂v terms in the momentum and in theenergy equations: v∼l , (15) ∂x (cid:12) (cid:12) dv ∇p 1 ν ∼l2(cid:12)(cid:12) ∂v(cid:12)(cid:12) , (16) = − +f + ∇·τ ∂x dt ρ ρ (cid:12) (cid:12) or(cid:12), mo(cid:12)re generally, for a betterstatistical evaluation, Navier-Stokesmomentum equation(10) (cid:12) (cid:12) ν2 ∂v ∂v 2 ∂v ∂v 2 ∼ x + x + y + y + ddt ǫ+ 12v2 =−ρ1∇·[(pv−v·τ)+c∇(ρǫ)]+f ·v l4 (cid:18) ∂y ∂z (cid:19)2 (cid:16) ∂x ∂z (cid:17) ∂v ∂v (cid:16) (cid:17) z + z . (17) Navier-Stokesenergy equation,(11) ∂x ∂y (cid:18) (cid:19) where the viscous stress tensor τ and the thermal con- Here, the problem relies in the evaluation of l. Being ductivity 1/ρ∇·[c∇(ρǫ)] terms are explicitly added, to be h the computational spatial resolution, and being L the solved together with thecontinuityequation,thekinematic scale length of the entire computational domain, h 6 l 6 equationandtheEoS.Itisimportanttonotethatthether- L. Because of the lack of any geometric information, the mal flux term includes two contribution: the first contribu- only physical scale lengths we know are those relative to tiondependsonthethermalgradient(∇ǫ)thatiscurrently the hydrostatic equilibrium (in the presence of an external usedforsolidsorforincompressiblefluids,whilethesecond force field): dp/ρf, as well as p/|∇p|, ρ/|∇ρ|, |v|/∇·v, one dependson thedensity gradient (∇ρ) correlated to the |v|/|∇×v|, etc.. R mass diffusivity, here not included in the continuity equa- Since,in3Dthenaturaltendencyisthedevelopmentof tion. The matrix element (α,β) of theviscous stress tensor smaller structures in a direct cascade process (Kolmogorov 1941a,b), some authors (Trampedach & Stein 2011) calcu- τα,β =ησα,β+ζ∇·v (12) late −1 and l= l−1 . (18) i σα,β = ∂∂xvα + ∂∂xvβ − 32δα,β∇·v. (13) Xi ! β α wherel referstovariousscalelengthsas(∂lnρ/∂r)−1, i η andζ arethedynamicfirst(shear)andsecond(bulk) (∂lnV/∂r)−1, (∂lnp/∂r)−1, (∇·v/|v|)−1, etc.. physical viscosity coefficients. Instead,noinformation wehaveabout v,sinceweonly In the present study,we simply consider ζ =0 and eq. knowv and c . flow s (7) as EoS. By definition, the physical kinematic viscosity In the viscous accretion disc modelling, the Shakura coefficient is ν =η/ρ. and Sunyaev Shakura (1972); Shakura& Sunyaev (1973) 4 G. Lanzafame parametrization ofturbulentviscosityislargely adopted.In Withoutanyconsiderationtotheexistenceortotheinvolve- this approach, thekinematic viscosity coefficient is mentofinternalenergylevelsforanidealflow,physicaldis- sipation transfers macroscopic ordered kinetic energy flows 1 ν = lv, (19) intoheat,thatisinmicroscopicchaotickineticenergyflows. 3 This means that elastic scattering collisional cross sections where both l and v are unknown. Assuming the flow foranidealgasarethoserightforareformulationofν.Fora isotropy,aKepleriantangentialkinematics,andthevertical gasmixture,themeanvalueoftheelasticscatteringimpact hydrostatic equilibrium, cross section l=αlH, (20) κ= Xiκi (23) where H ≃ rcs/vKepl, the local disc thickness, is the Xi local shortest macroscopic scale length, and αl 6 1 is a should be considered, where Xi =ni/ ini is the rel- scaling quantity.Without any isotropy assumption, αl >1. ativenumerical abundanceof thechemical species i. At thesame time, However,weneedtotakeintoaccounttPhetotalnumber ofmicroscopicmolecules-atomswithinthemacroscopiccross v=α c , (21) v s sectiondeterminedbythemixinglengthasineq.(18)which whereα 61.Wheneverv>c ,shockswoulddissipate arithmetically gives as a result a mixing length close to the v s the energy,reducing the velocity to subsonic. Hence, in the smallest scale length in the summation. For this reason, we Shakuraand Sunyaevapproach, prefer to computel as ν = 1αα c H =α c H, (22) l=min(l1,l2,l3,...,ln), (24) 3 l v s SS s the various l referring to p/|∇p|, ρ/|∇ρ|, |v|/∇ · v, with α <1 to befound. i SS |v|/|∇×v|, and so on. Pringle (1981) found 0.01 6 α 6 0.03 for 0.1 6 SS Hence,thetotalnumberofbarionscontainedwithinthe L/LE 61(LE =Eddingtonluminosity)asalowerlimitfor 3D mixing mass ρl3 is ρl3/µm , where m is the proton active galactic nuclei (AGN). For numerical simulations of H H mass and µis themean molecular weight. AGN,α ≈10−4−10−3isoftenadopted(Lanzafame et al. SS Since we need the molecular collisional cross section, 1998, 2008). Values for αSS ≈ 10−2 have also been found in 3D it is necessary to compute (ρl3/µm )2/3, that is the for theobserved protostellar objects (Hartman et al. 1998), numberof molecules within the volume l3H, powered to 2/3. whileforafitofFUOrionisobservedoutburst(Clarke et al. This number has to be multiplied with κ to get a statisti- 1990; Bell & Lin 1994; Lodato & Clarke 2004), α ≈ SS callyeffectivecollisional surfacecomposedofamultitudeof 0.001−0.003. Infullyionized discsindwarfnovae,thebest microscopic cross sections: observedevidencesuggestsatypicalα ∼0.1−0.4, whilst SS therelevantnumericalMHDsimulationsevaluateαSS ≈10 ρl3 2/3 times smaller. King et al. (2007) attribute this discrepancy κ (25) µm H to incorrect magnetic and boundary layer shortcomings in (cid:18) (cid:19) the computations. Nevertheless, this is not the conclusion To calculate ν, we need to divide this arithmetic term of the full story because Lanzafame (2008, 2009) showed byalength λ.Thislength decreases wheneverthelocal nu- that a well bound viscous accretion disc structure mod- merical density increases. Therefore λ∼n−1/3. ellingstronglydependsonseveralconditions:thekinematic As far as the velocity contribution of ν is concerned, ofthemasstransfer,γ,α andsoon.Forisothermalorfor we exclude not only v , but also the thermal ǫ1/2, being SS flow a quasi-isothermal thermodynamics, a disc is structurally this last strictly linked to the chaotic microscopic kinemat- bound evenfor α =0 because numerical dissipation only ics. Kinematic velocities of extraneous bodies - as in the SS is enough to produce a disc in shocks events, a result also original Von K´arma´n formulation - moving in the fluid are discussed by Sawada et al. (1987); Spruit et al. (1987) and, not considered. This exclusion is a consequence of the fact more recently,by Lanzafame (2010a) in its physical sense. thatthekinematicviscosity coefficient, mutuallycorrelated to the thermal conductivity and the diffusivity coefficients, are all expressions of the intrinsic physical property of the fluid.Hence, toconclude, our 4 ν AND C COEFFICIENTS AND MOLECULAR CHARACTERISTICS ρl3 2/3 ρl2 ν ≃ξ= κn1/3c = κc , (26) Different formulations of macroscopic physical dissipation µmH s µmH s (cid:18) (cid:19) rarely converge with each other and often include an arbi- being n=ρ/µm . trary parameter, to be evaluated. In addition, any corre- H Its2Dcounterpart,accordingtothesamealgebraiclog- lation to microscopic (molecular, atomic, nuclear) physical ical stepsis: properties is absent. Tothispurpose,weproposethefollowing evaluation of Σl2 1/2 Σl the physical kinematic viscosity coefficient ν, to be used to ν ≃ξ= µm κn1/2cs = µm κcs, (27) H H determine the viscous stress tensor τ in the Navier-Stokes (cid:18) (cid:19) equations. being Σ the2D mass density. Microscopic molecules, atoms, nuclei, haveknown elas- In these expressions, both the molecular/atomic µm H tic scattering impact cross section κ, useful to compute ν. and κ are included, as well as other macroscopic physical Molecular dissipation and macroscopic physical viscosity 5 quantities like Σ or ρ. This means that, the spatial com- c =3RT/2=3K T/m .Experimentally (c/ν)(ǫµ/c )≈ V B H V ponent of ν could strictly be not shorter than l and longer 1.3−2.5, instead of 1. The discrepancy is largely explained than h, according to other physical variables, especially Σ becauseofthefactthattheoreticallycisevaluatedconsider- or ρ. This formulation is free of any arbitrary parameter. ingauniformmolecularvelocitydistributioninsteadoflocal Moreover, in (26) and (27), the quantities (ρ/µm )lκ and molecular kinematic differences related to thepresence of a H (Σ/µm )κ are pure numbers both <1, otherwise the den- temperaturespatialgradient.Asaconsequence,considering H sity of the fluid is comparable or greater than the atomic ν =ξ, for an ideal gas, density. It is important to note that eqs. 26 and 27 are not c c K T exactlyequivalent.Ineq.26ν ∝l2 initsspatialdependence c=2ν V ≃2ξ V =3ξ B (29) ǫµ ǫµ ǫµm (as Prandtl’s eqs. 16, 17), while in eq. 27 ν ∝l in the same H spatialcontribution(asShakuraandSunyaev’seq.22).Eqs. so that, 26 and 27 should be considered strictly correlated either to 2 2 l κ c l c a 3D or to a 2D modelling, respectively. Indeed,densities ρ c≃2 ρ c V =3K T ρκ s (30) µ m s ǫ B µm ǫ andΣareformallycorrelatedbytheequivalenceoftheirnu- (cid:18) (cid:19) H (cid:18) H(cid:19) mericaldensitiesas(Σ/µmH)3=(ρ/µmH)2.SothatΣ≡ρl in 3D, and would bea very specific case. 2 2 Molecular or atomic collisional cross sections are as- c≃2 1 Σl κ c cV =3K T 1 Σlκcs (31) sumed, for the sake of simplicity, circular, without any dis- µ m s ǫ B µm ǫ H H (cid:18) (cid:19) (cid:18) (cid:19) tortion. Molecular dimensions are determined by the so in 2D. Since coefficients ν and c are comparable for called Van derWaals mean radius, definingthelimit where dilute gases, then theviscous and the thermal conductivity the microscopic molecular force potential becomes attrac- timescalesarealsocomparable.Thismeansthattheinertia tive, deviating from that of a non interacting free particle of matter tends toward a uniform kinematic and thermal relative to an ideal gas. In the case of fully ionized gas, the configuration with thesame time scales. repulsive Coulomb elastic scattering among head on collid- In the viscous computational fluid dynamics, currently ing ions determines theshortest classical impact parameter rpasrp ≃(4πǫ◦)−12Z1Z2e2/3KBT,relatedtoκasκ=πrp2. νcoarrnedlacteadrefronmoteoanclhyoatrhbeirt,raasriiltycpoarrraecmtleytrsihzeodu,ldbbute.aIlnsospniotet Z1eandZ2earethetwoeffectiveelectricchargesofthetwo ofthefactthatresultsinisothermalorinnearlyisothermal colliding ions, K is the Boltzman constant, T is the tem- B conditions could still be significant, the lack of any correla- perature and ǫ◦ is thedielectric constant of thevacuum. tion between ν and c is free of any physical meaning. Noticethatthedecreaseofthemixinglengthlinamore Intherestofthepaper,wemainlypayattentiontothe effectivecharacteristicscalelengthasaresultofthepresence physical viscosity. However, since now onwards, any con- of a multitude of small scale constituents, as shown in eqs. clusion referring to the role of the physical viscosity will (26,27),doesnotalterthemeaningofthekinematicviscos- also refer to the role of the thermal conductivity in a close ity coefficient role in the tensorial expression of the stress cause-effectcorrelation, whereanyanysignificantlocalspa- viscous tensor (eqs. 12, 13) in the Navier-Stokes equations tialderivativeinvolvingarelativemotionamongcontiguous describing fluid flows. Therefore, although we introduce lo- flowpartswillinvolveaabrakingandaviscousheating;any calphysicalpropertiesofthefluid,andalthoughturbulence local heating will involve larger pressure and temperature isnotafeatureoffluidsbutoffluidflows,thereformulation gradients and consequently a transport of energy and mass of the scalar ν within τ stays meaningful. α,β with comparable time scales; a heat transfer will involve a decrease of temperature and pressure spatial gradients. 4.1 Kinematic viscosity and thermal conductivity coefficients for dilute gases 5 VISCOSITY TESTS The kinematic viscosity coefficient ν and the thermal con- ductivity coefficient c are dimensionally identical. Both are Flowtransportanddampingofturbulencearetheonlytwo characterized by a scale length multiplied by a scale veloc- fields where the role of physical dissipation in the Navier- ity. Hence, both are transport coefficients. A relevant dif- Stokes equations is better emphasized. Therefore, in this ference between viscosity and thermal conductivity is that section,appropriatetestsonthephysicalviscosityefficiency while viscosity is activated only whenever a relative mo- willbedone.Inparticular,regarda2DSPHshocklessmod- tion occurs among contiguous flow elements, thermal con- ellingontheradialflowviscoustransportinanannulusring ductivityis an energy transport mechanism always existing anda2DSPHmodellingofBurger’sturbulence.γ =5/3will wheneveratemperaturegradientoccurseveninsteadystate be assumed throughout the tests. This low compressibility conditions. However, the two coefficients are always related regime is the less advantageous condition for viscous forces to each other because both explains the tendency of the againstpressureforces,aswellasagainstnumerical-artificial thermodynamic system toward a homogeneity and isotropy damping (Molteni et al. 1991; Murray 1996; Okazaki et al. status, smoothing out kinematic (ν) and thermal (c) local 2002),explainedbyamodest contributionof thebulkcom- spatial discrepancies. The ratio c/ν (R eif (1965)) is: ponent∇·vineqs.(12,13),notcompensatedbythehighcs, c c onν.Insuchtests,resultsonν =ξ,alsoincludingathermal = V , (28) conductivityc∝ν willalsobecomparedwiththoserelative ν ǫµ to where c is the molar specific heat of the gas at V constant volume which, for an ideal (monoatomic) gas is ν =csh (32) 6 G. Lanzafame ν =c l (33) s which are the two simplest analytical expressions for ν, where the characteristic length could be either at its minimum geometric value h, or coincident with the mix- ing length l. Despite their full mathematical meaning, for- mulations (32) and (33), like (eqs. 16, 17, 22), lack of any physical sense without any correlation to what matter the flow is made of. 0 5.1 2D radial viscous transport in a shockless isothermal annulus ring Theory on 2D shockless radial transport in a Keplerian an- nulus ring in a gravitational potential well (Pringle 1981) predicts that, from the Green function, the solution of the initial Keplerian mass distribution at time t=0 for Σ is: Σ(r,t=0)=mδ(r−r◦)/2πr◦ (34) for a ring having mass m and an initial radius r◦. The 0 0 0.5 1 1.5 02 0.5 1 1.5 2 solution,attimet,intermsofdimensionlessradiusx=r/r◦ r r and viscous time θ=12νtr−2 is ◦ Σ(x,t)=(m/πr−2)θ−1x−1/4exp[−(1+x2)/θ]I (2x/θ),(35) Figure 1.Examples offourcalculated massradialdistributions ◦ 1/4 asafunctionoftheviscoustimeθ.TwoαShakuraandSunyaev where I1/4 is the modified Bessel function. The action radialdistributionsareshown,aswellasaν=csh,withh=0.09, of viscosity is to spread out the entire annulus ring toward comparedwithourν=ξ radialprofile. a disc structure transporting most of the low angular mo- mentum mass toward the centre of the potential well and Prandtl formulation (15), do not yield any realistic radial transporting a smaller fraction of high angular momentum profile of mass distribution because of the too large viscos- mass toward theempty externalspace. ity, due to the very negligible spatial derivatives in the l Forpracticalcomputationalpurposes,sinceitisimpos- calculation, causing therapid accretion of theentire disc. sible to reproduce a delta Dirac function at time t = 0, it In spite of the low spatial resolution adopted for prac- is necessary to start numerical calculations from an initial tical purposes, results on radial shockless viscous transport massdistributionrelativetoasmallθvalue.Inourexample, herereported clearly showthat,without anyotherphysical we choose θ=0.017, a value comparable with that used by alternative,theShakuraandSunyaevformulationlookslike Speith & Kley (2003) for a radial viscous transport similar an appropriate expression for ν for 2D shockless disc struc- test. tures, since it uses a characteristic length H much shorter Fig. 1 shows a comparison among several radial distri- thanlalonga3rddimensionthatdoesnotexistin2D.How- butionsofthe2DmassdensityΣ(r,θ)forvariousθ.Because ever,physicalformulations as in eq. (26) or (27), where the oftheuniversalityofgraphicalrepresentationbyusingθ in- effective scale length (Σl/µm )κ is much shorter than l in stead of t, computational radial profiles at an assigned θ H thecaseofdiffusematter,couldbeavalidphysicalalterna- value have to correspond with each other, fitting the ana- tive,free of any arbitrary parameter. lytical solution (eq. 35). In particular, fig. 1 shows in the samepictureΣ(r,θ)profilesbothfortwoShakuraandSun- yaev kinematic viscosity ν =α c H, and for the ν =c h, SS s s 5.2 Damping of 2D Burger’s turbulent flow where h = 0.09, and for the newest ν = ξ, where our for- mulation for ν is used. This clearly shows that the viscous Statistical studies of turbulence normally involve hypothe- ν = ξ model correctly replicates the right radial shockless sesabouthomogeneityandisotropyonthedistributionand transport mechanism, without any local distortion. In this kinematicsofthespatialflow(Kolmogorov1941a,b).2Dtur- example, we used throughout the models, M◦ = 2·1033 g, bulence is relevant to understand large scale flows (Frish R◦ =1011cm,astypicalastrophysicalvalues.Weconsidered 1995; Kellay & Goldburg 2002). 2D turbulence schemati- aninitialdensityΣ◦ oftheorderof10−10 gcm−2,aninitial cally discusses either a ”forced steady state turbulence”, or sound velocity cs =5·10−2v◦, where v◦ =2π(GM◦/R◦)1/2 a”decayingturbulence”,ifanexplicitforcingtermisadded cms−1 isthenormalizationvalueforthevelocity,andagas in the momentum equation or not. Eddies of different sizes composed of pure molecular hydrogen. The thermal energy showing density and potential fluctuations in the flow nor- perunitmassǫiskeptconstant,beingdǫ/dt=0throughout mally characterize turbulence. theentiresimulation.Variationsoftheseparametersdonot 2Dphenomenologyissomewhatmorecomplexthan3D produceanydifference in theradial densitydistribution for phenomenology, even though computationally more conve- each θ, being the viscous time θ an absolute reference time nient (Tabeling 2002), being not derivable from simple di- for Σ. The ν =c h profile is shown as a further test on the mensionalarguments. 3Dturbulenceinvolvesscales smaller s radial profile, where eq. (32) has been considered in the ν than the trigger one and it is supposed to be hosted to a calculation. Instead, models adopting eq. (33), or even the direct enstrophy cascade from large to small scales (from Molecular dissipation and macroscopic physical viscosity 7 Figure2.(X,Y)plotsofdensitymapatvarioustimesT forthe Figure3.(X,Y)plotsasthoseofFig.2forthephysicallyviscous physicallynonviscousmodelwithh=0.05.64greytoneareused. modelν=csh. vX,vY tomogramsarealsoreported,showingvelocityfluctuation both duringthe initialturbulent phase andduringthe following dampingsubsequent phase. Turbulence cascade theory substantially predicts that the kinetic energy is contained in turbulent eddies, while thecascadeprocessforenstrophyissomewhatambivalentin small to large wavenumbers), where the mean kinetic en- 2D,accordingtoboundaryandinitialconditions.Thelarger ergy is transferred and mean enstrophy is conserved. In- theeddy,thelargeritskineticenergycontent,accordingtoa stead, 2D turbulence involves larger scales than the trigger power scaling law of eddy dimension, respecting theenergy one and it is supposed to be hosted to an inverse enstro- conservation law (kinetic + thermal) within the whole sys- phy cascade from small to large scales (from large to small tem.In2D,thesocalled”dualcascade”processdetermines wavenumbers).An inversecascade of energy and a contem- the formation of larger turbulent eddies up to the limit of porary direct cascade of enstrophy in 2D is called a ”dual” theentire spatial domain. cascade (Manz et al. 2009). Even though a physical viscos- The numerical experiment here carried out for 2D ity is usually considered, as responsible of flow damping of Burger’s turbulence studies the temporal evolution of the theNavier-Stokesequations,thephysicalcounterpartofar- damping of a chaotic gas inside a L = 5 size 2D squared tificialviscosityhasbeendiscussedinLeVeque(1992,2002); box, whose density and whose kinetic velocity are initially Fletcher (1991); Hirsch (1997). locally random, as an example of decaying turbulence. v X 8 G. Lanzafame Figure4.(X,Y)plotsasthoseofFig.2forthephysicallyviscous Figure5.(X,Y)plotsasthoseofFig.2forthephysicallyviscous modelν∼l2 ∂v . modelν=csl. ∂x (cid:12) (cid:12) (cid:12) (cid:12) andv valuesarewithintherange−5·10−4 and5·10−4 at turbulence dominates, the final configuration is character- Y thebeginningwhileh=0.05.Hencev ∼0andv ∼0sta- izedbyafluidstatistically atrest,wherefluctuationsinthe X Y tistically at time T =0. Throughout the models the initial velocity field, in the thermal energy and in the density are adimensional thermal energy perunit mass ǫ=1.86·10−7. reduced after some time toward a more general uniformity. The initial settings consists of an uniform thermal en- Figg.2to6showXY distributionofdensityatvarious ergy together with a statistical macroscopic homogeneous times T for various physical kinematic turbulent viscosities and isotropic spatial distributions on density and on veloc- ν,aswellasparallelvelocitytomogramsshowingv andv X Y ity components v , v . Being the macroscopic scale length atthesametime,fromtimeT =0toT ∼100.64greytones X Y much larger than thespatial resolving power, clear fluctua- for the density maps are shown between the minimum and tionexistsinthehomogeneityonthescaleresolutionlength. themaximum Σvaluesforeach plot.That is, each plot has This is enough to ignite the initial turbulent kinematics, itsownminimumandmaximumvalueforΣ.Thenormaliza- withlargervelocity,thermalenergyanddensityfluctuations tionvalueforthedensityisΣ◦ =10−10gcm−2anddensities as larger are pressure forces from the beginning. Damping areinitiallycomputedmultiplyingΣ◦ timesarandomnum- viscous effects always work since the beginning. Hence, as berbetween0and1.Wedistinguishtwophysicalregimes.In a decaying test, after the initial turbulent condition, where the first one, from the beginning to T ≈25, 2D turbulence Molecular dissipation and macroscopic physical viscosity 9 Figure6.(X,Y)plotsasthoseofFig.2forthephysicallyviscous Figure7.(X,Y)plotsasthoseofFig.2forthephysicallyviscous modelν=ξ foraH2 gas.Σ◦=10−10 gcm−2. modelν=ξ foraO2 gas.Σ◦=1.6·10−9 gcm−2. dominates because small eddies grow up. Linear size of ed- spatialresolutionlengthispossible,ofcourse,butnotprac- dies,initiallyoftheorderofh,growsuptoseveralhvalues. tical. The reduction of thetime step, together with a much Therelevanceoftheartificialdissipation-ν ≈(0.1−0.4)c h larger computer memory needed involve very long compu- s -isevidentinthephysicallynonviscouscase,becauseofthe tational time, even for 2D simulations. We adopted a low low spatial resolution adopted.From T ∼25 onwards,clear spatial resolution, since h/L=10−2. This involves that re- differences both on the density spatial distribution and on sultsinherenttoallmodelsarequiteviscous.However,these the velocity component fluctuation appear, distinguishing results, even though conditioned by the intrinsic numerical the efficacy of the various adopted physical viscosity coef- damping effect, clearly show how ξ correctly works, free of ficients, since the physical viscosity better works increasing any dependenceon h, or on any arbitrary parameter. c , as well as enlarging the eddies if ν ∝ l or ∝ l2. Plots Toshowtheroleofotherchemicalspecies,throughtheir s clearly show that ν = c l is the most dissipative model in meanmolecularweight,wealsoshowinFig.7resultsregard- s so far as Σ◦ =10−10 g cm−2. In this example, the physical ingthedampingof2DturbulenceforagasofmolecularO2 damping relative to ν =ξ is comparable with that relative forΣ◦ =1.6·10−9 gcm−2attimeT =0.Theseresultshave to thePrandtl’s ν ∼l2|∂v∂x|. tobecomparedwiththoseforagasofmolecularH2withthe A reduction of theintrinsic damping by decreasing the same initial numerical density (Fig. 6), whose Σ◦ = 10−10 10 G. Lanzafame also because of a more relevant heat thermal conductivity suppressing local thermal differences caused by the viscous heating. At the same time, the viscous heating is so strong thattheconsequentincreaseofpressureforcestimebytime activates again the velocity fluctuations in the flow after T ∼ 50. Hence, a different evolution of velocity fluctuation isdeveloped,dependingon theinitial conditions forv, ǫ,ν, Σ.Thestudyofthe2Dviscousdampingofthe2Dkinemat- ics of a so densefluid does not need to concern us. 6 3D ACCRETION DISC IN A CLOSE BINARY We compare results for 3D stationary disc structures both in high compressibility (γ = 1.01) and in low compressibil- ity (γ = 5/3) with the aim of getting a physically well- bound accretion disc around the compact primary star in a LMCB. However, these comparisons are thought to show which physically viscous disc structure, in theShakura and Sunyaev formulation (eq. 22), in the Prandtl formulation (eqs. 16, 17), and in the more simple formulations (eqs. 32 and 33), better compare with that relative to ν = ξ. The adopted spatial resolution length is h = 5·10−3 through- outthesimulations, beingthemutualseparation ofthetwo stars normalized to 1. Thecharacteristicsofthebinarysystemaredetermined by the masses of the primary compact white dwarf and of itscompanion starandtheirseparation.Wechosetomodel asysteminwhichthemassM1 oftheprimarycompactstar and the mass M2 of the secondary subgiant star are both equalto1M⊙ andtheirmutualseparationisd12 =106 Km. The injection gas velocity at L1 is fixed to v ≃ 130 Km inj s−1 while the injection gas temperature at L1 is fixed to T◦ = 104 K, taking into account, as a first approxima- tion, the radiative heating of the secondary surface due to radiation coming from the disc. Supersonic kinematic conditions at L1 are discussed in Lanzafame (2008, 2009); Lanzafame et al.(2000,2001),especiallywhenactivephases ofCB’sareconsidered.Thereferenceframeiscentredonthe primary compact star, whose rotational period, normalized to2π,coincideswiththeorbitalperiodofthebinarysystem, beingthevelocitiesnormalizedtov◦ =[G(M1+M2/d12)]1/2. Results of this paper are to be considered a useful test to Figure8.(X,Y)plotsasthoseofFig.2forthephysicallyviscous check whether disc structures (viscous and non) show the modelν=ξ foraH2 gas.Σ◦=10−8 gcm−2. expected behaviour. We simulated the physical conditions at the inner and at theouter edges as follows: g cm−2 at T = 0. Results for O2 are influenced by the gas a) inneredge: molecular crosssection κ,insofarasthenumericaldensity thefreeinflowconditionisrealizedbyzeroinggasflowinside ρ/µm isrelevantintheξcalculationtogetaνcomparable a sphere of radius 10−2, centred on the primary compact H orhigherthananyartificialoranynumericaldissipation.In star.Althoughdiscstructureanddynamicsarealterednear this example, being the numerical densities for the molecu- theinneredge,thesealterationsarerelativelysmallbecause lar O2 at the beginning the same as those of Fig. 6 for H2, they are balanced by a high gas concentration close to the ν =ξ is so effective that viscous dissipation quickly damp- inneredge in supersonic injection models. ensbothanyinitialvelocityfluctuationandflowprogression b) outer edge: toward any spatial homogeneity. gas flow from L1 towards theinterior of theprimary Roche AnevenstrongerviscousflowdampingisshowninFig. Lobe is simulated by constant gas pressure, density, and 8foradensergasofmolecularH2,whoseΣ◦ =10−8 gcm−2 thermal energy per unit mass, as well as a constant veloc- atT =0.Viscousmoleculardamping,whereν =ξ,appears ity in a small conic region having L1 as a vertex and an more effective for higher density flow modelling. As an ex- apertureof ∼57◦.Theradiallength ofthissmall volumeis ample, in Fig. 8 viscous dissipation is so relevant that any ∼ 10h. The initial injection particle velocity is radial with localflowmotionisstronglysuppressedsincethebeginning, respect to L1. Local density at the inner Lagrangian point