Discussion of Transport Coefficients

The dotted black lines in Figs. 8.1 show that the two-parameter fit according to Eq. (8.1) with Eqs. (8.6) is very good for all coefficients of restitution α and all global area fractions φ0

presented. In a large region in the middle of each system the fit is hardly distinguishable from the directly measured heat flux. For higher inelasticities (α . 0.7−0.8) and lower densities (φ0 . 0.1−0.2) the fit becomes worse. This is not surprising because in these cases the gradients increase significantly leading to a failure of the gradient expansion, Eq. (8.1).

In the left graph of Figs. 8.2 we show both transport coefficients κ (full red lines) andµ/ε (full dark green lines) as a function of the coefficient of restitutionαfor a very wide range of coefficients of restitution0.2 ≤α≤0.999for moderately dense systems,L_x = 20,L_y = 25, N = 256, φ₀ = 0.4. The left graph of Figs. 8.3 shows the same observables for0.5 ≤ α ≤ 0.999for a total area fractionsφ₀ = 0.2. The blue lines showκ_F from a fit to Fourier’s law.

The dashed orange and light green lines show the globalκ_JR0 andµ_JR0/εfrom Eqs. (8.2).

Both transport coefficients,κandµ/ε, extracted from our data are non-monotonic inα. Start-ing at inelastic systems and increasStart-ing α, the thermal conductivity κ first decreases, goes through a minimum and then sharply increases again to tend to a non-zero constant in the elastic limit for all densities studied. The other transport coefficient is divided byεbecause we assumeµ∝ε. Its behavior depends on the density. For moderate densities (Figs. 8.2), coming from inelastic systems, µ/εfirst increases with increasing α, reaches a maximum and finally decreases again. This maximum is in the vicinity of the minimum of κbut clearly not at the same value ofα. For lower density systems (Figs. 8.3),µ/εfirst decreases when the inelastic-ity is reduced, reaches a minimum, and sharply increases again. It is not clear whether it will reach a maximum and decrease again or keep increasing and possibly reach a finite value.

The insets show a magnification of the quasi-elastic region and reveal that the measured coeffi-cient of thermal conductivityκis in good agreement with the predictions by Jenkins and Rich-man [JR85] for α ≥ 0.99and moderately dense systems (Figs. 8.2). For higher inelasticities

0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Global transport coefficients κ and µ for Lx=20, Ly=25, N=256 κµ/ε

Figure 8.2.: The left graph shows the transport coefficients κandµ/εfor systems withLx = 20, Ly= 25, andN = 256, corresponding to a line densityλ= 10.24and a total area fractionφ0 = 0.4.

We compare the values obtained from our simulations with the corresponding predictions from Ref.

[JR85], Eqs. (8.2), for0.2≤α≤0.999. Note, however, that the theory in Ref. [JR85] assumes quasi-elastic systems – therefore agreement is not expected for values ofαthat are not close to one. The inset of the left graph shows a magnification of of the quasi-elastic regime. The right graph quantifies the importance ofµby showingR1andR2as defined in Eqs. (8.7) and (8.8), and discussed in the text.

significant deviations are observed. Since Jenkins and Richman [JR85] assume quasi-elastic systems we cannot expect their results to be valid for higher inelasticities. Still, their theory yields the correct order of magnitude of the measured κ for allα. As to µ_JR0 the deviations toµare quite strong. The discrepancy ofµ(more than one order of magnitude in most cases) is unclear to us. For lower densities (Figs. 8.3) it has not been possible to obtain the trans-port coefficients for sufficiently elastic systems to judge whether or not the measured data will converge to the predictions by Jenkins and Richman [JR85].

Soto et al. [SMR99] have performed molecular dynamics (MD) simulations to study µfor a very dilute granular gas on a vibrated plane under gravity. They made use of the non-monotonic behavior of the temperature within such systems to extract µat that position where the tem-perature gradient vanishes – a method that works in the presence of gravity only, even though their results are virtually independent of the magnitude of gravity. This suggests that it should be valid for zero-gravity environments, too. In comparisons we find that their coefficientµis always much smaller than ours, up to one order of magnitude. This suggests that the evaluation of µis a subtle task andµmight depend on global system parameters – similarly to the local equation of state in chapter 7. Other reasons might be that zero-gravity environments behave qualitatively different than gravitational ones and, probably most importantly, the simulations of Ref. [SMR99] apply to very low densities only and not to0.1≤φ₀ ≤0.4. Furthermore, the difficulties in our fitting procedure increase considerably asα → 1. In this limit the temper-ature gradient and the density gradient both become linear inxand thus proportional to each other, see Figs. 6.12 and 6.13 (page 51) and Eqs. (6.4) and (6.5) on page 57. Consequently, it is not possible forα →1to determinetwoparameters from the fit unambiguously. Finally, we

0.5 0.6 0.7 0.8 0.9 1

α

10^-4 10^-3 10^-2 10^-1

Global transport coefficients κ and µ for Lx=40, Ly=25, N=256

κµ/ε κ_F κ_JR µ_JR0 0/ε

0.97 0.98 0.99 1

0.03 0.06 0.09 0.12

0.5 0.6 0.7 0.8 0.9 1

α

0 5 10 15 20

|χF|/|χ|-1 vs. α for Lx=40, Ly=25, N=256

0.5 0.6 0.7 0.8 0.9 1

0.03 0.1 0.3

|µ∆|/|κ_FδT|

Figure 8.3.: Same as Fig. 8.2 but for systems withLx= 40, corresponding to a total area fraction of φ0= 0.2. Compared to Fig. 8.2 the dip inκis shifted towards elastic systems andµ/εdoes not decrease again. Of course, the plainµstill vanishes forα→1. The right graphs are qualitatively similar to the ones in Fig. 8.2.

expect the transport coefficients to be x-dependent and not constant as assumed in our fitting procedure.

To estimate the degree of inhomogeneity of the transport coefficients inxwe have divided the sample box into an inner and an outer part and fitted the heat current using data from either half of the box only. The scattering of the data from the inner and outer part is shown in Fig.

8.4 (indigo and magenta lines, respectively) and provides a rough measure for the effects of inhomogeneous transport. The results for the outer half (magenta lines) are almost identical to the ones for the full system (black lines). This is also true for the inner part (indigo lines) except for the weakly inelastic systems for which the absolute value of the heat flux in the middle of the sample is so small that statistical fluctuations dominate. This shows that the determination of the transport coefficients clearly needs further investigation, especially in the quasi-elastic regime.

Let us now turn to the importance of the new transport coefficientµ. The right graphs of Figs.

8.2 and 8.3 show the quantitative measures (R₁andR₂) of the contribution to the heat flux that isnot proportional to the temperature gradient. The main graphs showR₁ as defined in Eq.

(8.7) and the insets showR₂as defined in Eq. (8.8). From the main graphs as well as from the insets we see that the contribution that is not proportional to the temperature gradient becomes negligible in the elastic limit, as expected. The insets further show that this contribution is never greater than 30% for all systems studied and that it is more important for larger inelasticities.

The peaks inR1(main graphs) approximately indicate the border line up to which Eq. (8.1) is a good approximation to the measured heat flux. As already mentioned, fits to Fourier’s law are hardly distinguishable from the ones to Eq. (8.1) for quasi-elastic systems. For moderately inelastic systems, shown to the right of the peaks in the right graphs of Figs. 8.2 and 8.3, the

0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 α

10^-3 10^-2 10^-1

Global transport coefficients κ and µ for Lx=20, Ly=25, N=256

κ, 0<x/L_x<0.33 µ/ε, 0<x/L_x<0.33 κ, 0<x/L_x<0.165 µ/ε, 0<x/L_x<0.165 κ, 0.165<x/L_x<0.33 µ/ε, 0.165<x/L_x<0.33 κ_JR

0 from Ref. [JR85]

0.9 0.92 0.94 0.96 0.98 1

0.10 0.15 0.20

µ(α)/ε κ(α)

Figure 8.4.:Transport coefficientsκandµ/εfor the same systems as in Fig. 8.2 (Lx= 20,Ly = 25, andN = 256, corresponding to a line densityλ = 10.24 and a total area fractionφ0 = 0.4). In addition we show here the values for the transport coefficients obtained for the inner (indigo) or outer part (magenta) only. The inset comparesκfrom our simulations to the value of Ref. [JR85].

two-parameter fit works well while Fourier’s law already fails – leading to a high value ofR₁. For even lower inelasticities, shown to the left of the peaks, Eq. (8.1) starts to fail, too, such that neither description is able to faithfully reproduce the directly measured heat flux. This causes R₁to decrease.

Distributions in Vibrated Systems

Finally, we examine the local velocity distributions of the driven granular gas in the stationary state. We distinguish between the local distribution

f_x(x, v_x) := 1 ρ(x)

Z

R

dv_y f_stat(x, v_x, v_y) (9.1) at positionxof the velocity componentv_xin the direction of the driving and the distribution

f_y(x, v_y) := 1 ρ(x)

Z

R

dv_xf_stat(x, v_x, v_y) (9.2) of the velocity componentv_y perpendicular to the direction of the driving. By definition, these velocity distributions are normalized to unity.

In order to determinef_xandf_yfrom the simulation we use two different methods. The first one extracts them directly according to their definitions (9.1) and (9.2) fromfstatin the following way: We count the number of particlesf(r, v_x, v_y, t)|Vr|dv_xdv_yat timetin cellVrwith an x-component of the velocity betweenv_xandv_x+dv_x andy-component of the velocity between v_y andv_y+dv_yand average over a long time interval. One disadvantage of this method is that it yields a smeared-out velocity distribution, which is spatially averaged over the width of the stripVxcentered aroundx. This is of practically no importance in the middle of the simulation box, but strongly disturbing for resolving the subtleties which occur close to the driving walls and are presented in Sec. 9.1. The second way of measuring the velocity distributions avoids this problem: To getf_i (i = xor y) we keep track of all those particles which pass the line parallel to the y-axis at position xwithin a very long time interval of lengthτ and whoseith component of the velocity lies in a small interval of width∆v around v_i. Let us enumerate these particles byniand denote their velocity components byvx⁽ⁿⁱ⁾andvy⁽ⁿⁱ⁾. Then one has (in the limitsτ → ∞and∆v→0)

ρ(x)f_i(x, v_i) = 1 τ∆vLy

X

n_i

1

|v⁽ⁿ_x ⁱ⁾|, (9.3)

because fori=x, resp.i=y, the right-hand side of (9.3) is equal to Z

R

dvy |F_x(x, v_x, v_y)|

|v_x| resp.

Z

R

dvx |F_x(x, v_x, v_y)|

|v_x| . (9.4)

75

Here,F_x(x, v_x, v_y) :=v_xf_stat(x, v_x, v_y)denotes thex-component of the (differential) current density in the stationary state at positionxof particles with velocity componentsv_xandv_y. As compared to the first method of measuring velocity distributions, this one is also statistically more effective for determining rare events, such as the high-velocity tails in Sec. 9.2. However, as far asf_x(x, v_x)is concerned for smallv_x, the second method is inferior to the first, because (9.3) assigns a large weight to the relevant events and therefore amplifies statistical fluctuations, too.

9.1. Effects of the Discontinuity at a Driving Wall

Particle-number conservation at a driving wall requires the incoming particle flux at the wall to be equal to the outgoing flux. For the velocity distribution f_x this implies the boundary condition [BRMM00]

f_x(∓L_x/2, v_x) =Θ(±v_x−v_dr)

1∓v_dr v_x

f_x(∓L_x/2,−v_x+v_dr), (9.5) which must hold for allvx >0. HereΘ(z) := 1forz≥0, respectivelyΘ(z) := 0forz <0, denotes Heaviside’s unit-step function and v_dr = 1 in our units chosen. (Note that we use a driving velocity v_dr = v₀ here instead ofv_dr = εv₀ as used in chapter 6). The boundary condition (9.5) relates the distribution f_x of velocities prior to a collision with the wall to the one after a collision with the wall. In contrast, the velocity distributionf_ymust obey the usual reflection symmetry inv_y everywhere in the system, that is

f_y(x, v_y) =f_y(x,−v_y) (9.6) for all|x| ≤L_xand allv_y >0.

We measured the velocity distribution fx(x, vx) at 25 different positions in a moderately in-elastic system with coefficient of restitution α = 0.9. Fig. 9.1 shows the rescaled velocity distribution f˜_x(x, v_x/p

T_x(x)) := f_x(x, v_x)p

T_x(x), measured using method one. At the driving walls f˜_x is seen to obey the boundary condition (9.5). When moving from a driving wall towards the center of the system, the gap inf˜_xgets gradually smeared out, andf˜_xbecomes more and more symmetric. Clearly, even for this moderately inelastic system f˜_xdoes not scale for differentx. The two extreme cases forx =−Lx/2andx = 0are shown together in Fig.

9.2. The data forf˜_x at the wall was obtained with the second method for measuring velocity distributions. This result is in agreement with Direct-Simulation-Monte-Carlo results in Ref.

[BRMM00] and Molecular-Dynamics simulations in Ref. [BT02b].

-3 -2 -1 0 1 2 3 -0.5 -0.4

-0.3 -0.2

-0.1 0

0.1 0.2

0.3 0.4

0.5

0 0.2 0.4 0.6 0.8 1 1.2

PSfrag replacements

x/Lx

vx/p Tx(x)

Figure 9.1.:Rescaled velocity distributionsf˜xfor differentx. [System parameters:α= 0.9,Lx= 20, λ= 10.24(N = 256,φ0= 0.4)]

-3 -2 -1 0 1 2 3

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

left wall

middle

PSfrag replacements

v^x/p T^x(x)

Figure 9.2.: Two rescaled velocity distributionsf˜xfrom Fig. 9.1: at the left wall and in the middle of the system.

-3 -2 -1 0 1 2 3 0

0.1 0.2 0.3 0.4

f_x in the middle f_y at the wall f_y in the middle Gaussian

-1 -0.5 0 0.5 1

0.25 0.3 0.35 0.4

f_x middle f_y wall f_y middle Gaussian

~

~

~

~

~

~

PSfrag replacements

vi/p Ti(x) f˜x

f˜y fx˜ fy˜

Figure 9.3.: Combined plot of rescaled velocity distributions of the same system as in Figs. 9.1 and 9.2: f˜xin the middle of the sample (green crosses),f˜y at a driving wall (red squares), andf˜y in the middle (blue triangles). Also shown is a centered Gaussian with unit variance (black solid line). The inset shows the central part of the main graph with the (1433/2964/1340) data points being smoothed (running average over (41/85/38) data points).

9.2. Non-Universal High-Velocity Tails

Our main result for the velocity distributions is the non-scaling and multiformity of their tails.

In order to observe these phenomena, extensive simulations for capturing rare events are re-quired and the data have to be analyzed on a logarithmic scale. In contrast, on alinearscale the rescaled velocity distributions seem to collapse approximately, as was observed previously by e.g. [GZBN97, MSS04]. As an example we show in Fig. 9.3 the rescaled distributions f˜_y(x, v_y/p

T_y(x)) := f_y(x, v_y)p

T_y(x)in the middle of the sample (x = 0) and at a driving wall (x=±Lx/2) for the same moderately inelastic system (α = 0.9) as in Figs. 9.1 and 9.2.

Clearly, the curves are reflection-symmetric inv_y. For comparison, we have also included the symmetric distribution f˜_xin the middle of the system. Deviations are visible for small veloc-ities only, as the inset of Fig. 9.3 shows. The corresponding Gaussian (solid line) in Fig. 9.3 suggests that the velocity distributions are close to but not identical to a Maxwellian. The ap-proximate data collapse, observed on this level of accuracy in Fig. 9.3, even continues to hold if the coefficient of restitution of the gas is varied in a not too large extent. In considerably more inelastic systems, such as forα = 0.5, this is not true any longer. For example, the peak off˜_xmeasured in the middle of the system would be considerably broader and flatter than the ones off˜_y in the center and at the wall (not shown).

-10 -5 0 5 10

Figure 9.4.: Rescaled velocity distributions f˜x in the middle of the sample (green crosses), f˜y at a driving wall (red squares), andf˜y in the middle (blue triangles) as in Fig. 9.3, but here smoothed data are shown on a semi-logarithmic scale and velocities of much higher absolute values are included. The solid line is a centered Gaussian with unit variance. Different parts of the figure represent different systems with decreasing global area fractionφ0(left to right) and decreasing coefficient of restitutionα (top to bottom). The remaining system data areLx= 20,λ= 10.24(N = 256) for the left column and Lx= 50,λ= 9.6(N = 240) for the right column. The insets show|ln ˜fi|⁻¹on a double-logarithmic scale to determine the decay exponentβfrom (9.7). The dashed line is the best linear fit tof˜xin the middle of the sample, the solid line tof˜yin the middle, and the dotted line tof˜y at the wall. Part a) corresponds to the system of Fig. 9.3.

In contrast, Fig. 9.4 a) shows the same data of Fig. 9.3 on asemi-logarithmicscale and includes also velocities of much higher absolute values. From this figure it is evident that scaling does not hold in the high-velocity tails of the distributions either. Similar observations were made before in e.g. [BT02b, vZM04, BRM03]. The type of decay in the high-velocity tails is differ-ent forf˜xandf˜y, and also depends on the position in the sample, the coefficient of restitution α, and the global area fractionφ₀. This is illustrated by examples of different systems in Fig.

9.4. In the insets we show |ln ˜f_i|⁻¹ on a double-logarithmic scale in order to determine the decay exponentβdefined by

ln ˜f_y 0, v_y/√

T_y(0)|vy|→∞

∼ − v_y/√

T_y(0)

^β (9.7)

forf˜_yin the middle of the sample. The exponent is defined accordingly forf˜_yat the wall andf˜_x in the middle of the sample. For the moderately inelastic systems withα= 0.9in the first row of Fig. 9.4, the asymptotics has been clearly reached. We note thatβis different for the different distributions, and also depends on the global area fractionφ₀. Upon loweringα(top to bottom in Fig. 9.4) and/or decreasing the global area fractionφ₀ (left to right), the tails off˜_yget more and more populated, that is,β decreases. In very simple terms, this may be understood from the fact that (i) the largest typical velocities are always of the order ofv_dr = 1(see also Fig.

9.5) and (ii) thatv_dr/p

T_i(x)increases up to 20 with decreasingαand decreasingφ₀. Hence, a Maxwellian velocity distribution would not be able to supply enough probability to particles with velocities of the order of v_dr, instead higher-populated tails are needed. This argument suggests different behavior in different velocity regions so that the distributions cannot be fitted to the functional form (9.7) over the entire range of velocities [vZM04]. Indeed, such a behavior can be seen in Fig. 9.4 d), e) and f). The final asymptotics could not always be deduced from the simulations, even though our data include velocities which are up to 40 times bigger than the appropriate granular velocities √

T_i. This applies tof˜_y in the middle of the sample in part e), where we suspect that the final asymptotics has not been reached. An asymptotic analysis of f˜_xin the middle of the sample is even more problematic due to particles that reach the system center from a driving wall without undergoing a collision. These particles give rise to the side peaks off˜_xin parts e) and f), which have prevented us from determiningβin these cases.

Fig. 9.4 contains smoothed data as in the inset of Fig. 9.3, and the scales of the horizontal axes are determined by the square root of the appropriate granular temperatures. For comparison, Fig. 9.5 shows the unsmoothed data corresponding to Fig. 9.4 f), plotted directly versusv_i (in units ofv_dr). The side peaks ofp

T_x(0)f_x(0, v_x)in the middle of the sample (green crosses in Fig. 9.5) are due to the above-mentioned particles which fly from a driving wall to the system center without undergoing a collision. Concerning the distribution of they-component of the velocity in the middle of the sample (blue triangles in Fig. 9.5) we find for the vast majority of particles that |v_y|is less than 80% of the driving velocity. It is in the region of this highest velocity observed for p

T_y(0)f_y(0, v_y) in the middle of the sample (blue triangles), where pT_x(0)f_x(0, v_x)in the middle of the sample (green crosses) changes abruptly its slope.

-2 -1 0 1 2 10^-6

10^-5 10^-4 10^-3 10^-2 10^-1 10⁰

PSfrag replacements

v

_i

pTx(0)fx(0, vx)

√Ty(0)fy(0, vy)

√Ty(Lx/2)fy(Lx/2, vy)

Figure 9.5.:Velocity distributions for the system in Fig. 9.4 f), but with unsmoothed data and without rescaling of the horizontal axis.

The exponentβhas also been determined experimentally in a strongly driven gas so that gravity effects are small [RM00]. A value β ≈ 1.55 was measured for a gas with coefficient of restitutionα≈0.93. It was found to be remarkably independent of the global area fractionφ₀, which was varied from0.05to0.3. We have also simulated the system of [RM00] with zero gravity (not shown) and reproduced their value forβ. Even though the experimental data cover a wide range of velocities, some of the interesting phenomena discussed in Fig. 9.4 cannot be observed in this range. For the same reason the significance of the particular valueβ ≈ 1.55 should not be overestimated.

Simulations of a driven granular gas in a circularly shaped box also yield stretched Gaus-sian tails [vZM04]. In addition, evidence is given that β depends only on the coefficient of restitution and the average ratio of the number of particle-wall collisions and particle-particle collisions. For homogeneously driven systems the tail of the velocity distribution was the-oretically predicted [vNE98] to be governed by the decay exponent β = 3/2. Simulations of homogeneously driven systems [MSS01, BT03] observed the exponent β = 3/2only for unrealistically low values of the coefficient of restitution.

Spheres with Coulomb Friction

Incorporating Coulomb friction turns out to be quite challenging. For this reason we restrict our analysis of frictional particles to homogeneously driven systems. The driving mechanism has been specified in Sec. 2.2, see page 12. Binary collisions with Coulomb friction are modeled using the Walton model, cf. Sec. 2.1.2, page 10.

In a homogeneously driven system, the long-time averages of all hydrodynamic fields are ho-mogeneous. Instead it is interesting to study the mean translational and rotational temperatures, their stationary ratio, and their relaxation to the steady state. In this chapter we will present and discuss an analytic mean field calculation to obtain differential equations for the full time evolution of the translational and rotational temperature. Technical details of the calculation can be found in appendix C.

Since the analytic calculations are very elaborate we will introduce and probe various sim-plifications to the collision rules for Coulomb friction (the Walton model, Eq. (2.8)). These simplifications (labeled models B through D) will be introduced in Sec. 10.2. Model A is the

Since the analytic calculations are very elaborate we will introduce and probe various sim-plifications to the collision rules for Coulomb friction (the Walton model, Eq. (2.8)). These simplifications (labeled models B through D) will be introduced in Sec. 10.2. Model A is the

Im Dokument Stationary Properties of Driven Granular Gases (Seite 81-0)