Visco-rotational shear instability of Keplerian granular flows

TL;DR

This paper investigates the linear stability of viscous Keplerian granular flows using a rheological model, revealing an instability driven by shear rheology and differential rotation that impacts planetary rings and protoplanetary disks.

Contribution

It introduces a new rheological instability mechanism in dense granular flows under Keplerian rotation, highlighting conditions for its onset.

Findings

Instability occurs when viscosity grows faster than the square of shear rate.

The instability influences the dynamics of planetary rings.

It may affect granular flows in protoplanetary disks.

Abstract

The linear stability of viscous Keplerian flow around a gravitating center is studied using the rheological granular fluid model. The linear rheological instability triggered by the interplay of the shear rheology and Keplerian differential rotation of incompressible dense granular fluids is found. Instability sets in granular fluids, where the viscosity parameter grows faster than the square of the local shear rate (strain rate) at constant pressure. Found instability can play a crucial role in the dynamics of dense planetary rings and granular flows in protoplanetary disks.

Equations

$$\rho \left\{\frac{\partial }{\partial t}+V_k\frac{\partial }{\partial x_k}\right\}{V}_i=-\frac{\partial P}{\partial x_i}+\rho \frac{\partial \Phi }{\partial x_i}+\frac{\partial {\tau }_{ik}}{\partial x_k},$$

$${\tau }_{ik}=\eta {\dot{\gamma }}_{ik},{\dot{\gamma }}_{ik}=\partial V_k/\partial x_i+\partial V_i/\partial x_k,$$

$$\partial V_k/\partial x_k=0.$$

$$\eta =\eta (P,\xi ),\xi =\sqrt{{\dot{\gamma }}_{ik}{\dot{\gamma }}_{ik}/2}.$$

$$r{\Omega }^2=-\frac{\partial \Phi }{\partial r},$$

$$\left(r\frac{{\partial }^2\Omega }{\partial r^2}+3\frac{\partial \Omega }{\partial r}\right)\bar{\eta }+r\frac{\partial \Omega }{\partial r}\frac{\partial \bar{\eta }}{\partial r}=0,$$

$$\Phi (r,z)=\frac{GM}{(r^2+z^2{)}^{1/2}}$$

$$\Omega (r)={\Omega }_0{\left(\frac{r}{r_0}\right)}^{-q},{\Omega }_0={\left(\frac{GM}{r_0^3}\right)}^{1/2}.$$

$$\frac{\partial \ln \bar{\eta }}{\partial \ln r}=q-2.$$

$$\mathbf{V}=\bar{\mathbf{V}}+{\mathbf{V}}^{\prime },P=\bar{P}+P^{\prime },\eta =\bar{\eta }+{\eta }^{\prime }.$$

$${\bar{V}}_{\varphi }(r)=r_0{\Omega }_0+{\frac{\partial (r\Omega )}{\partial r}|}_{r_0}(r-r_0)+...$$

$$x=r-r_0,y=r_0(\varphi -{\Omega }_{0}t),$$

$$A=\frac{r_0}{2}{\frac{\partial \Omega }{\partial r}|}_{r_0},B=-{\Omega }_0-A,$$

$${\bar{V}}_y(x)=2Ax$$

$$\frac{\mathrm{D}}{\mathrm{D}t}{V}_x^{\prime }-2{\Omega }_0{V}_y^{\prime }$$

$$\frac{\mathrm{D}}{\mathrm{D}t}{V}_y^{\prime }-2B{V}_x^{\prime }$$

$$\frac{\mathrm{D}}{\mathrm{D}t}{V}_z^{\prime }$$

$$G_P\equiv {\left(\frac{\partial \eta }{\partial P}\right)}_{\xi },G_S\equiv \frac{1}{\rho }{\left(\frac{\partial \eta }{\partial \xi }\right)}_P.$$

$$\frac{{\eta }^{\prime }}{\rho }=G_P\frac{P^{\prime }}{\rho }+G_S\left(\frac{\partial V_y^{\prime }}{\partial x}+\frac{\partial V_x^{\prime }}{\partial y}\right).$$

$$\left(\begin{array}[]{c}{\bf V}^{\prime}({\bf r},t)\\ P^{\prime}({\bf r},t)/\rho\\ \eta^{\prime}({\bf r},t)/\rho\end{array}\right)\propto\left(\begin{array}[]{c}{\bf u}({\bf k},t)\\ -{\rm i}p({\bf k},t)\\ -{\rm i}\mu({\bf k},t)\end{array}\right)\exp\left({\rm i}{\bf r}{\bf k}(t)\right)~{},$$

$${\dot{u}}_x(t)$$

$${\dot{u}}_y(t)$$

$${\dot{u}}_z(t)$$

$$0$$

$$\mu (t)$$

$$\omega =\pm {\left({\bar{\kappa }}^2-W^2\right)}^{1/2}+\mathrm{i}\left(W-\nu k^2\right),$$

$${\bar{\kappa }}^2=\left(-4B\Omega -4A^2{G}_{\gamma }{k}_x{k}_y\right)\frac{k_z^2}{k^2-4A{G}_P{k}_x{k}_y},$$

$${\sigma }_A$$

$${\sigma }_P$$

$${\sigma }_S$$

$${\bar{\kappa }}^2>W^2,$$

$${\bar{\kappa }}^2<W^2,$$

$$G_P<0.$$

$$G_S>0.$$

$${\left(\frac{\partial \ln \eta }{\partial \ln \xi }\right)}_P>2.$$

Full text

Disk of solid particles rotating around central gravitating object is an important class of granular flows widely occurring in nature. Among those are planetary and exoplanetary rings, debris disks around young stars, or even areas of protoplanetary disks where dust particles accumulate and form dense granular material. These flows, occurring at different scales, often have several common features: solid particles rotate on nearly Keplerian orbits, highly inelastic particle collisions can easily dissipate kinetic energy, and self-gravity of granular material can be neglected in comparison to the gravitational potential of the central object. Granular flows normally collapse into thin disks, where particle number density increases and in some cases the flow can be described using a fluid model with “granular viscosity”.

It is known from the accretion disk theory that differentially rotating viscous flows can be unstable Lightman and Eardley (1974); Shakura and Sunyaev (1976). Indeed, it has been shown that viscous instability sets in when the increase of surface density leads to the decrease of the local viscosity Lin and Bodehnheimer (1981); Ward (1981); Lukkari (1981). In this case, smallest density bump leads to the enhanced angular momentum transfer and corresponding accretion process. Hence, mass accumulation at the outer edge of the perturbation leads to further increase of density. Viscous instability can operate in optically thick disks, where the viscous stress is proportional to radiation pressure. However, phenomenological tests reveal somewhat uncommon character of the instability, which even when occurred, provides insignificant growth rates for linear perturbations.

The second alternative energy source in Keplerian granular flows is the viscous overstability Kato (1978); Blumenthal et al. (1984) that is thought to be a primary mechanism for the development of some of the observed structures in dense planetary rings. This axisymmetric pulsational instability occurs in granular flows, where the derivative of kinematic viscosity with respect to the surface density is positive and exceeds some critical value Ward (1981); Papaloizou & Lin (1988); Schmid and Tscharnuter (1995). Thus, compressible epicyclic response leads to viscous overcompensation and growth of density-spiral waves due to increase of the viscous stress in the compressed phase. Later viscous overstability has attracted considerable interest including its non-axisymmetric Borderies et al. (1985); Longaretti & Rappaport (1985); Spahn et al. (2000); Schmidt et al. (2001); Latter & Ogilvie (2006) as well as nonlinear saturation properties Schmit & Tscharnuter (1999); Schmidt & Salo (2003); Latter & Ogilvie (2010).

The key to the investigation of granular flows around gravitating objects is a proper account for the particle collision effects. Kinetic description of particle collisions has been successful in modelling properties of rapid and dilute granular flows. Still, kinetic approach may fail due to the scale separation problem between granular and flow time-scales and inelasticity of particle collisions. In fact, it is known that a detailed theoretical description of granular flows should deal with number of specific features: granular gases are intrinsically non-equilibrium systems with non-Maxwellian distribution functions that in some cases can reveal non-local of even non-Markovian character (see Ref. Goldhirsch (2003) and references therein). Still, granular flows can be studied using hydrodynamic equations that can describe collective phenomena including different types of instabilities, thermal convection Khain and Meerson (2003); Pontuale et al (2016), behaviour of granular gas mixtures, or clustering Olafsen and Urbach (1998); Brey et al (1998); Goldhirsch (1999); Serero et al. (2006).

Significant advances in the understanding of the dense granular fluids have been made recently. It seems that a wide range of dense granular flows can be unified into a rheological model that permits formulation of a local constitutive equation Jop et al. (2006); Forterre and Pouliquen (2008); Boyer et al. (2006). In this local rheological model granular phenomenology is employed to define how fluid viscosity depends on pressure, as well as strain tensor of the flow. Thus, granular flow can be described by incompressible non-Newtonian fluid model, where strain tensor is solely due to the velocity shear of the flow. We employ this model for the description of astrophysical flows, where individual dust granules can be highly porous particles colliding with a low restitution parameters. In this limit dense granular flow can exhibits “fluid” properties even at moderate values of particle volume fraction.

In the present paper, we study the linear stability of viscous Keplerian flow around a gravitating center, taking into account rheological aspects of the viscous stress tensor. Our incompressible model includes pressure and shear rheology since they both affect linear stability of spiral waves. We identify unstable axisymmetric modes analytically and analyze non-axisymmetric instability numerically.

Physical model. The dynamics of an incompressible viscous flow rotating around central gravitating object can be described by the Navier-Stokes equation:

[TABLE]

where $\rho$, $P$ and $V_{i}$ are density, pressure and velocity of the flow, respectively. We neglect self-gravity and assume that $\Phi$ is the gravitational potential of the central object. The viscous stress tensor $\tau_{ik}$ can be calculated using the strain rate tensor

[TABLE]

in incompressible limit it is reduced to a shear strain tensor:

[TABLE]

To describe the dissipative properties of the dense granular flow we employ rheological fluid description implying the existence of a local constitutive equation. Indeed, it has been shown recently, that granular fluids can be described using the specific form of the non-Newtonian fluids (see Ref. Jop et al. (2006) and references therein). In this limit viscosity of granular fluid $\eta$ depends on both, pressure as well as the second invariant of the strain rate tensor $\xi$:

[TABLE]

This frictional visco-plastic constitutive law has been tested successfully in laboratory experiments and is thought to be a general model describing dense granular flows in “fluid” regime Forterre and Pouliquen (2008). The “fluid” regime of dense granular flows in laboratory is realized for a narrow range of granular volume fraction, defined as the ratio of the volume occupied by the grains to the total volume. Still, the granular rheology used here may also work for lower density systems where the coefficient of restitution is low.

Alternative interpretation of the rheological model set by Eq. (4) can be obtained within the assumption of microscopic turbulence. Indeed, Boussinesq eddy viscosity hypothesis assumes that turbulent viscosity parameter can be calculated using the strain rate tensor (see Eq. 2). In such limit, eddy viscosity can vary due to the variation of the intensity of microscopic turbulence, depending on the pressure or local velocity shear of the flow.

Steady state. Let us consider axisymmetric stationary differentially rotating viscous flow in the cylindrical coordinates with constant pressure $\bar{P}$ and density $\bar{\rho}$. Azimuthal velocity of the background depends on the angular velocity of the differential rotation $\bar{V}_{\phi}=r\Omega(r)$. The radial and azimuthal components of the Navier-Stokes equation of the stationary state in polar frame reads as:

[TABLE]

where

[TABLE]

is the gravitational potential of the central object with mass $M$. Assuming thin disk model ($z^{2}/r^{2}\ll 1$) we derive rotationally supported steady state where the gravitational potential of central object sets Keplerian profile of the angular velocity:

[TABLE]

Here $r_{0}$ is some fiducial radius used to parameterize the steady state and $q=3/2$. Hence, using Keplerian angular velocity into the Eq. (6) we can derive radial profile of the viscosity parameter in equilibrium:

[TABLE]

Interestingly, Rayleigh stability criterion in rotating fluids $\partial_{r}(r^{2}\Omega(r))>0$, or $q<2$, indicates that in steady state, viscosity parameter should be a decreasing function of radius: $\partial_{r}\bar{\eta}<0$. Hence, Eqs. (8) with radially homogeneous pressure and density form the globally stable granular Keplerian flow that can be used for the local linear stability analysis.

Local linear analysis. To study the linear dynamics of dense granular flows we split the velocity, pressure and viscosity parameter into the background and perturbation components:

[TABLE]

We employ local shearing sheet approximation, where the flow curvature effects can be neglected and the differential rotation is reduced to the plane shear flow Goldreich and Lynden-Bell (1964); Tevzadze et al. (2003, 2010). In this limit we expand azimuthal velocity

[TABLE]

and use local approximation to neglect higher order terms with respect to $(r-r_{0})/r_{0}$. Hence, introducing the local Cartesian frame co-rotating with the disk matter at the fiducial radius $r_{0}$

[TABLE]

and using standard form of the Oort constants

[TABLE]

we can calculate steady state velocity

[TABLE]

that describes the radial shear of the azimuthal velocity due to the differential rotation of the flow.

Hence, equation governing the linear dynamics of the perturbations in local shearing sheet frame can be reduced to the following:

[TABLE]

where $\nu=\bar{\eta}/\rho$, ${\rm D}/{\rm D}t\equiv{\partial/\partial t}+2Ax{\partial/\partial y}$ and $\Delta=\partial^{2}/\partial x^{2}+\partial^{2}/\partial y^{2}+\partial^{2}/\partial z^{2}$ and the radial gradient of viscosity parameter is neglected in the local approximation: $\partial\bar{\eta}/\partial x=0$.

To describe rheological properties of the flow we employ a general form of the local constitutive equation and introduce pressure $G_{P}$ and shear $G_{S}$ rheology parameters as follows:

[TABLE]

Assuming that the rheological parameters of the granular fluid can be considered to be locally constants we can calculate linear perturbation of the viscosity as follows:

[TABLE]

Introducing Fourier expansion of the spatial variables in shearing sheet frame

[TABLE]

where ${\bf k}(t)=(k_{x}(t),k_{y},k_{z})$ and $k_{x}(t)=k_{x}(0)-2Ak_{y}t$, we can derive the system of equations governing the linear dynamics of incompressible perturbations in time:

[TABLE]

where $\dot{\psi}(t)$ stands for the time derivative of the variable $\psi(t)$ and $k^{2}(t)=k_{x}^{2}(t)+k_{y}^{2}+k_{z}^{2}$. Equations (21) pose a complete initial value problem that can be solved numerically. However, to get more insight into the stability properties of the system we derive an approximate dispersion equation.

Stability analysis. Dispersion equation of the ODE system (21) can be derived in the case of rigid rotation ($A=0$). However, we employ adiabatic approximation when time dependent mode frequency can be introduced and linear perturbations can be expanded in time as: $\psi(t)\propto\exp(-{\rm i}\omega(t)t)$. In this limit we assume that frequency depends on time only through the shearing variation of wave numbers: $\omega(t)=\omega({\bf k}(t))$. Thus, the dispersion equation leads to:

[TABLE]

where $\bar{\kappa}$ sets epicyclic frequency in rheological flows:

[TABLE]

and $W=\sigma_{A}+\sigma_{P}+\sigma_{S}$ with

[TABLE]

Here $\sigma_{A}$ describes the shear flow transient amplification due to the differential rotation of the flow, while $\sigma_{P}$ and $\sigma_{S}$ describe the effects of pressure and shear rheology, respectively.

In the rigidly rotating Newtonian fluids ($G_{P}=G_{S}=0$) solution reduces to the classical spiral wave dumped by constant viscosity: $\omega=\pm 2\Omega_{0}|k_{z}/k|-{\rm i}\nu k^{2}$.

The existence of growing modes can be seen in the case of differentially rotating flows. Eq. (22) shows that the necessary condition for the growth of linear perturbations in differentially rotating granular fluids is $W>0$. Therewith, the character of the perturbation growth depends on whether rheological stress can destabilize epicyclic balance or not:

[TABLE]

Axisymmetric perturbations. Eq. (22) is rigorous in describing the stability of axisymmetric modes with $k_{y}=0$. In this limit transient amplification is absent ($\sigma_{A}=0$), and we can analyze rheological modifications of the spiral waves.

For the purpose of direct comparison with the viscous instabilities we neglect shear rheology ($G_{S}=0$) and analyze the effect of pressure rheology parameter. Then the necessary condition of the perturbation growth reduces to:

[TABLE]

This in turn indicates that the viscous overstability developing at $\partial\eta/\partial\rho>0$, i.e., $G_{P}>0$ is an intrinsically compressible mechanism that is absent in the incompressible limit.

In the opposite limit, when pressure rheology can be neglected ($G_{P}=0$), we recover new type of growth mechanism that originates from the shear rheology of the granular fluid:

[TABLE]

For better understanding we reformulate growth criteria as $\sigma_{S}=-AG_{S}k_{x}^{2}>\nu(k_{x}^{2}+k_{z}^{2})$. Hence, unstable modes are nearly uniform in the vertical direction $|k_{z}/k_{x}|\ll 1$. Using Eqs. (4,13,18) and local value of incompressible strain rate $\xi(r_{0})=-2A$ we may rewrite the shear rheology instability condition in a more general form:

[TABLE]

Thus, the shear rheology of the fluid leads to the visco-rotational instability when the granular viscosity parameter increases faster than the square of the shear (strain) rate.

In general, when pressure and shear rheology effects are comparable, necessary condition of instability can be reduced to the following: $\sigma_{P}+\sigma_{S}>\nu k^{2}$. Here we introduce the viscous cut-off wave-number $k_{\nu}$ that defines length-scales that normally dissipate during one rotation period: $\Omega_{0}=\nu k_{\nu}^{2}$. Hence, dynamically active modes are located in the $k/k_{\nu}<1$ area of the spatial spectrum.

The growth rates of linear axisymmetric perturbations are shown on Fig. 1. The growth mechanism due to pressure rheology favors large-scale perturbations ($k_{x}/k_{\nu}\ll 1$, panel A), while shear rheology instability operate at small radial scales ($k_{x}\sim k_{\nu}$, panel B). In all cases most unstable modes are nearly uniform in the vertical direction $k_{z}/k_{\nu}~{}\ll~{}1$. The growth rates of the visco-rotational instability set by the shear rheology are asymptotically higher at wave-numbers larger than the cut-off wave-number $k_{\nu}$. However, at length-scales shorter than the granular dissipative scales the very validity of the rheological model breaks down leading to the modification of the visco-rotational instability, a process that we do not address in the current paper.

Non-axisymmetric perturbations. Linear dynamics of non-axisymmetric modes can be analyzed through Eqs. (20-24), or numerical solution of the initial value problem (see Eqs. 21). Fig. 2 shows the growth rates in ($k_{x},k_{y}$) plane. Shearing sheet modes are drifting in this plane due to the background shear ($k_{x}=k_{x}(t)$). Thus, the non-axisymmetric modes have some finite time before reaching viscous scale $k_{\nu}$, where they are dumped due to a viscous dissipation. It seems that the pressure rheology parameter introduces leading-trailing asymmetry of the linear modes: leading modes grow higher for $G_{P}>0$, and trailing modes for $G_{P}<0$. Therewith, positive pressure rheology decreases the growth rates of the shear rheology instability, while the negative pressure rheology enhances it. Fig. 3 shows results of the numerical calculations of Eqs. (21). The energy of spiral waves is shown at different values of azimuthal wave-number. Figure illustrates the transient character of the growth of non-axisymmetric modes.

To get more insight into the nature of the instability we derive dynamical equation in the limiting case of vertically uniform perturbations ($k_{z}=0$) and shear rheology ($G_{P}=0$). We can reformulate Eqs. (21) for the horizontal velocity circulation:

[TABLE]

where ${\rm curl}({\bf u})_{z}=k_{x}(t)u_{y}-k_{y}u_{x}$ is the linear perturbation of the horizontal vorticity and Oort’s constant $A=-q\Omega_{0}/2$. Hence we may conclude that visco-rotational shear instability of horizontal vorticity perturbations occurs at $G_{S}>0$ in differentially rotating flows with angular velocity decreasing outwards, and at $G_{S}<0$ if $q<0$.

Summary. We present the new type of instability in a rheological viscous dense granular flows rotating around a central gravitating object. The incompressible visco-rotational instability originates from the shear rheology of the granular fluid. The instability operates on small scales and differs in principle from the known viscous instabilities due to the pressure rheology of viscous Keplerian flows. The mathematical formulation of the problem is set to demonstrate fundamental nature of the found instability. We adopt minimal model approach, showing that degrees of freedom necessary for this instability to develop are 3 dimensionality and supercritical shear rheology of the flow. The instability occurs in flows where the viscosity parameter has a positive steep gradient with respect to the local shear velocity. Unstable modes have small radial and large vertical scales, indicating the possibility of instabilities for narrow azimuthal rings (ribbons).

The visco-rotational shear instability can be simply described using the pressure-vorticity balance. For instance, anticyclonic vorticity perturbations to the Keplerian flow lead to local increase of the pressure. When this vorticity increase leads to the increase of the viscosity and corresponding accretion rate, pressure will increase even more, setting the linearly runaway process. A similar process will occur with cyclonic vorticity at pressure minima, for which a viscosity decrease will result in further the flow pressure decrease.

The visco-rotational shear instability may lead to a nonlinear saturation at higher amplitudes, or to the delocalization of the local constitutive relation and deveopment of non-local structures due to the specific properties of granular media Kamrin & Koval (2012). We speculate that the instability analyzed here can play a crucial role in the dynamics of dense planetary rings, as well as promote structure formation in protoplanetary disks in the areas of high dust to gas ratios.