Molecular dissipation in the nonlinear eddy viscosity in the Navier-Stokes equations: modelling of accretion discs

0.7 MB·English

by  Giuseppe Lanzafame

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