Investigation of a Vorticity-preserving Scheme for the Euler Equations
Darryl Seligman Karim Shariff



















































Darryl Seligman
darryl.seligman@yale.edu


Astronomy Department, Yale University,
    New Haven, CT 06511, USA


NASA Ames Research Center
Moffett Field, CA 94035, USA


















We investigate the vorticity-preserving properties of the compressible, second-order residual-based scheme, “RBV2”. The scheme has been extensively tested on hydrodynamical problems, and has been shown to exhibit remarkably accurate results on the propagation of inviscid compressible vortices, airfoil-vortex interactions on a curvilinear mesh, vortex mergers in an astrophysical accretion disk, and the establishment of a two-dimensional inverse cascade in high-resolution turbulent simulations. Here, we demonstrate that RBV2 sustains the analytic solution for a one-dimensional shear flow. We assess the fidelity by which the algorithm maintains a skewed shear flow, and present convergence tests to quantify the magnitude of the expected numerical dispersion. We propose an adjustment to the dissipation in the algorithm that retains the vorticity-preserving qualities, and accurately incorporates external body forces, and demonstrate that it indefinitely maintains a steady-state hydrostatic equilibrium between a generic acceleration and a density gradient. We present a novel
numerical assessment of vorticity preservation for discrete-wavenumber vortical modes up to the Nyquist wavenumber. We apply this assessment to RBV2 in order to quantify the extent to which the scheme preserves vorticity for the full Euler equations.  We find that RBV2 perfectly preserves vorticity for modes with symmetric wavenumbers, i.e., k_x=k_y, and that the error increases with asymmetry. We simulate the dynamical interaction of vortices in a protoplanetary disk to demonstrate the utility of the updated scheme for rendering astrophysical flows replete with vortices and turbulence.  We conclude that RBV2 is a competitive treatment for evolving vorticity-dominated astrophysical flows with minimal dissipation.

























§ INTRODUCTION
 
Simulations of high Reynolds number flows that respect the evolution of vorticity, and minimize its numerical diffusion, are of key importance for many applications of astrophysical fluid dynamics.  Numerical schemes that accurately treat vorticity are of special importance for high-value problems such as blade-vortex interactions in aeronautical engineering. Any numerical investigation of flows with high Reynolds numbers, where turbulence and vortices are expected to dominate the dynamics, inherently require a numerical method that accurately evolves the vorticity, as well as the density, momenta, and energy in the flow. In the astrophysical context, hydrodynamical or magnetohydrodynamical turbulence is thought to be the main driver of the dynamics in convective domains such as Jupiter's atmosphere <cit.>, in accretion disks <cit.> and in protostellar disks <cit.>, to name a few examples.


The Euler equations, and the vorticity equation derived from them, have no dissipative terms.  However, when solving the Euler equations numerically on a finite mesh, a certain amount of dissipation is required to represent the transfer of kinetic energy and enstrophy from  resolved to  unresolved scales.
On the other hand, numerical simulations generally suffer from excessive energy and enstrophy dissipation such that their spectra begin to rapidly fall-off at length scales substantially larger than the numerical grid size.


This is due to numerical dissipation, which is either inherent in the method (such as upwinding), or explicitly added to guarantee a stable solution. Schemes used to capture shock-waves employ
nonlinear dissipation. Particle-based schemes also tend to diffuse vorticity due to the dissipation associated with the adopted smoothing kernel. A classical approach to decrease diffusivity, and thereby increase the effective Reynolds number in the flow solution, is to increase the resolution of the simulation. Since hydrodynamical simulations are computationally expensive, it is of interest to investigate alternative approaches to controlling unphysical vorticity diffusion. 

Theoretically, since enstrophy cascades to small scales, long  simulations  must allow for unresolved vorticity structures. Therefore, simulations where turbulence and vortices are expected to dominate the dynamics should aim to model the effect of the unresolved scales on the resolved scales. The first model of this kind was proposed by  <cit.>, and from comparisons with direct numerical simulations it is known that it is too dissipative. Recently there has been an effort to develop less dissipative models <cit.>. Models  should ideally dissipate enstrophy  only at the scales near the grid cut-off where non-linear interactions between the resolved and unresolved motions are important. 

<cit.> introduced the concept of vorticity-preservation. A scheme is said to be vorticity-preserving if the evolution equation for discrete vorticity looks the same as the continuous vorticity equation discretized by central difference operators. This requires that the numerical operators, denoted by a tilde, for curl and gradient satisfy (curl) ( grad) = 0, which holds if the derivative operators for the different spatial directions commute.  It also requires that upwind and/or dissipation terms have zero curl. They demonstrated that the Lax-Wendroff-Ni scheme <cit.>, exactly preserves the discrete vorticity for the Euler equations in the case of linearized acoustics in a uniform medium.

<cit.> presented the numerical scheme RBV2, which extended the concept of vorticity preserving numerical simulations to the linearized (about a uniform flow) and full compressible Euler equations. The method incorporates additional dissipation into a central difference scheme, which by construction, respects the vorticity evolution.  The RBV2 algorithm has been extensively tested on hydrodynamical problems. <cit.> showed that the scheme is able to perfectly sustain the two-dimensional, compressible inviscid vortex proposed by <cit.>. <cit.> compared the same test problem with commonly used astrophysical codes, and found that RBV2 remarkably outperformed a Godunov-type solver using third-order Runge–Kutta timestepping and the weighted essentially nonoscillatory WENO3 <cit.> reconstruction scheme, a van-Leer advection based code and a spectral code. <cit.> demonstrated that the RBV2 scheme performs favorably on a linearly advected vortex and demonstrated its fidelity on the computation of airfoil-vortex interactions on a curvilinear mesh. <cit.> also demonstrated that, in the comparison codes, RBV2 was the only code among those tested that was capable of establishing and sustaining an inverse cascade in an unforced two-dimensional turbulence simulation. Motivated by the extremely efficient test results for RBV2, <cit.> implemented the RBV2 scheme in the context of astrophysical accretion disks. They found that it ran at a much higher effective Reynolds number than the Godunov method for comparable computational effort, based on a test of vortex merger with a background shear. They concluded that the algorithm may have important ramifications for simulations of astrophysical environments where turbulence was expected to dominate the dynamics, and where strong shocks are unlikely to be present.

It is of interest to note that there are other methods to solve the Euler equations that are capable of running at  high Reynolds numbers, and which minimze the diffusion of vorticity. <cit.> found that the spectral code Dedalus <cit.> was able to  sustain the inviscid vortex, although it took ten times longer to run than the RBV2 algorithm as a consequence of the very small time step required to run at very high Reynolds number.  

<cit.>, however, noted that when RBV2 is used to model flows that contain significant advection diagonally through the grid, spurious short-wavelength perturbations may manifest in
the crosswise direction. In order to dissipate these perturbations, while maintaining the desirable vorticity-preserving quality of RBV2, <cit.> implemented high-order spatial filters as proposed by <cit.>. The details of this implementation may be found in Section 2 of <cit.>.

In this paper, we continue to explore the potential of the RBV2 scheme, and specifically, we conduct an investigation of the vorticity-preserving qualities of the algorithm. We present several variations of a shear flow as a  hydrodynamical test, all ran to the same simulation time. These tests provide a linear wave convergence test, and quantify the magnitude of the numerical dispersion in the method. We propose an update to the algorithm to incorporate external forces, which is necessary for most astrophysical applications, such as numerical simulations of a protoplanetary disk or a planetary atmosphere. We implement this update, and demonstrate that the new algorithm is able to sustain a shear flow in the presence of a vertical hydrostatic equilibrium. We present  a novel assessment of vorticity preservation, which we apply to the algorithm. We conclude with  a simulation of the dynamical interaction of seven vortices in a protoplanetary disk, to demonstrate the utility of the updated scheme for rendering turbulent astrophysical flows. 

The plan of this paper is as follows. In  2, we present an overview of the algorithm, extending the formal analysis of its merits (and shortcomings), and proceeding from the work of <cit.>. In  3 we briefly present the numerical methods used to implement the scheme. In  4 we present the numerical and analytical results for RBV2 on a one-dimensional shear flow.  In  5, we present the results of numerical experiments on a skewed shear flow. In  6, we propose an update to the algorithm to incorporate external body forces and demonstrate that the updated scheme sustains a hydrostatic equilibrium. In  7, we provide a numerical quantification for the dissipation of vorticity for a set of vortical modes with wavenumbers varying from 1 to the Nyquist wavenumber. In 8, we perform shearing sheet simulations of vortex interactions in a protoplanetary disk as a demonstration of the updated scheme. In 9, we discuss the ramifications of these results, outline future prospects, and conclude.  



§ OVERVIEW OF THE ALGORITHM


Following the nomenclature of <cit.>, we write the two-dimensional Euler equations as,


    ∂_t  w +∂_x  f+∂_y  g=0   ,

for three general vectors w, f and g. The exact form of these three vectors depends on the equation set to be solved; <cit.> consider three sets of vectors for the cases of pure acoustics, the linearized Euler equations, and the full Euler equations. Since they will be used frequently in this paper, we provide the latter two here. For the Euler equations linearized around a uniform flow (u_0,v_0), the vectors are defined as:

    w   = [ p; u; v ]  ,
         f   = [ pu_0+u; uu_0+p;   vu_0 ]  ,
              g   = [ pv_0+v;   uv_0; vv_0+p ]  ,

whereas for the full Euler equations, they are given by, 

    w   = [   ρ; ρ u; ρ v; ρ E ]  ,
         f   = [        ρ u;  ρ u^2 + p;       ρ uv; (ρ E + p)u ]  ,
              g   = [       ρ v;      ρ uv;  ρ v^2 +p; (ρ E +p)v ]  .

For simplicity, we adopt the polytropic equation of state, p=c_s^2ρ^γ. For the remainder of this paper, we will further specialize to an isothermal equation of state, i.e., γ=1, with a uniform sound speed.


The RBV2 scheme may be cast in the form, 

    D_0 μ_0  w+ D_1 μ f+ D_2 μ g=h_1 D_1 (Φ_1 q̃)+h_2 D_2 (Φ_2 q̃)   ,

where μ_0 = μ^2 is a grid-averaging operator and μ will be defined below.
The terms on the right-hand-side of (<ref>) represent the added dissipation, h_1,h_2 being positive parameters tending towards zero with increased resolution. Φ_1 and Φ_2 are matrices with the same eigenvectors as the Jacobians A=df/dw and B=dg/dw respectively. The general difference operators are D_0 , D_1, and D_2, which represent the time, x and y numerical derivatives respectively, and the general averaging operators are μ and μ_0. q̃ is, in the general case, a vector involving discrete functions of w, f and g, although for the RBV2 algorithm, q̃=D_0μ w+D_1 f+D_2  g. For the remainder of the paper,  μ=μ_1μ_2, D_0=1/Δ tΔ̅, D_1=1/δ xδ_1μ_2 and D_2=1/δ yδ_2μ_1. These operators are defined on a Cartesian mesh with x_j=jδ x,  y_k=kδ y,  with time steps t^n=nΔ t, for an arbitrary primitive variable ψ, as,


    (δ_1 ψ)_j+1/2,k   =ψ_j+1,k-ψ_j,k
    
        (δ_2 ψ)_j,k+1/2   =ψ_j,k+1-ψ_j,k
    
        (μ_1 ψ)_j+1/2,k   =1/2(ψ_j+1,k+ψ_j,k) 
    
        (μ_2 ψ)_j,k+1/2   =1/2(ψ_j,k+1+ψ_j,k) 
     
        (Δ̅ψ)^n+1   =3/2ψ^n+1 - 2ψ^n +1/2ψ^n-1  .



<cit.> primarily consider Ω≡∇×ρu, as the quantity of interest and refer to it as the vorticity.  By vorticity-preservation they mean that the evolution equation for its discrete analog be identical to the equation obtained when the same difference operators are applied to the continuous evolution equation for Ω.  We adopt the same terminology and refer to ω⃗≡∇×u as the proper vorticity.

The discrete vorticity evolution equation  comes from taking the curl of the momentum part of Equation <ref>,


    D_0μ_0 curl_w+D_1μcurl_f+D_2 μcurl_g
     = h_1D_1curl_1+h_2D_2curl_2  ,


where


    curl_w   =D_1  w^(3) -D_2  w^(2)
    curl_f   =D_1  f^(3) -D_2  f^(2)
    curl_g   =D_1  g^(3) -D_2  g^(2)
    curl_m   =D_1(Φ_m q̃)^(3)-D_2(Φ_m q̃)^(2) , m = 1, 2 .



Theorems 1–5 in <cit.> prove that the scheme given by Equation <ref>, with the condition,

    h_1Φ_1   =h A 
     
        h_2 Φ_2   =hB   ,

for some constant h, is vorticity preserving for the case of pure acoustics and the  Euler equations linearized around a uniform advecting flow. Without repeating the details of each proof, we note that in general, they adhere to the following structure: (i) Define the scheme, (ii) Solve for the vorticity evolution (Equation <ref>) within the assumptions and (iii) Show that the all terms in the right hand side (RHS henceforth) of Equation <ref> reduces to zero. The extra dissipation here, by construction, demands that the discrete solution to the Euler equations respects the  vorticity equation to machine precision.



Theorem 6 in <cit.> states that the Scheme <ref> with Condition <ref>, preserves vorticity for the linearized Euler Equations (Equation <ref>), if the discrete operator,


    Λ = μ - h_1D_1Φ_1-h_2D_2Φ_2,

is non-singular. For the linearized system, the discrete vorticity evolution, Equation <ref>, reduces to 


    (D_0μ_0+u_0D_1μ+v_0D_2 μ) ω̃=h_1D_1curl_1+h_2D_2curl_2  ,

where the numerical (proper) vorticity ω̃≡ D_1v-D_2u.
They show that the Scheme <ref> is equivalent to the matrix equation,


    Λq̃ = 0   .

Therefore, if Λ is invertible, this yields q̃=0. They conclude that, since  the RHS of Equation <ref> is then zero,  the scheme is vorticity-preserving. However, we note that if q̃=0,  the RHS of the original scheme (<ref>) is identically zero; in other words, the additional dissipation in the original scheme vanishes.  This raises two questions: (i) why  is it necessary to introduce the additional dissipation in the first place, and (ii) would its absence mean that shocks would not be captured without spurious oscillations?


Theorem 7 in <cit.> extends Theorem 6 to the case of the full Euler equations, and leads to the same issue of the added dissipation being zero.
It was this issue that led us to numerically investigate the algorithm; see <ref>.
There we demonstrate that RBV2 does not preserve vorticity in general, which implies that the added dissipation terms are non-zero in general.

However, there is now a large amount of numerical evidence that the RBV2 algorithm is considerably less dissipative than many schemes.  We therefore conducted the following investigations to elucidate the extent to which the scheme preserves vorticity. 



§ SOLUTION OF THE IMPLICIT SYSTEM

In this section, we describe the numerical method by which we solve the implicit system (<ref>).  For a more detailed description, we refer the reader to <cit.>.  Equation (<ref>) is written as a residual equal to zero, i.e., R( w) = 0.
We solve this equation by stepping the following equation to steady-state (in fictitious time) using the technique proposed by <cit.>:


    ∂ w/∂ t^* + R( w)=0   ,

where t^* is a fictitious time. The residual is split into its convective and dissipative parts (Q and D, respectively),


    R( w)= Q( w)+ D( w)  .

The solution vector at subsequent time steps w^n+1 is found by iterating over the following algorithm, where the superscript k symbolizes successive iterations for w,


    w^n+1,0= w^n   ,


and successive iterations are


    w^n+1,k= w^n-α_kΔ t^*[ Q^k-1+ D^k-1]  .


The initial convective and dissipative residuals are calculated from the previous time step, Q^0= Q[ w^n] and D^0= D[ w^n]. The kth residuals are calculated using Q^k= Q[ w^n+1,k] and D^k=β_k D[ w^n+1,k]+(1-β_k) D^k-1. The coefficients α_k and β_k are chosen to maximize the stability interval along the imaginary and negative real axis of each Fourier mode of the solution. <cit.> reported  that  effective choices for α and β are α_1=1/3, α_2=4/15,  α_3=5/9, α_4=1 and β_1=1, β_2=1/2, β_3=0, and β_4=0. For stability, in all simulations presented in this paper, we use a global timestep that is one third of that given by the Courant condition, e.g. that Δ t = 1/3Δ x/√(u^2+v^2+c_s^2), where u and v are the maximum velocities at any point in the grid, and we use three sub-iterations of the scheme at each step. 





§ ONE-DIMENSIONAL SHEAR FLOW

In this section we examine both analytically and numerically, the fidelity with which RBV2 applied to the full Euler equations (<ref>) maintains the one-dimensional shear flow:


    u(x,y)   =Usin(my) 
     
        v(x,y)   =0 
    ρ(x,y)   =ρ
    
        p(x,y)   =c_s^2ρ  .


Since the equation of state is isothermal, we may neglect the internal energy equation. The resulting Jacobian matrices are,


    A= [         u         1         0; u^2+c_s^2        2u         0;        uv         v         u ]  , 
         B= [         v         0         1;        uv         v         u; v^2+c_s^2         0        2v ]  .


We demonstrate that the dissipation term in the scheme given by Equation <ref>, and therefore its discrete curl, is identically zero for the one-dimensional shear flow.  The only nonzero part of q̃ involves the y derivative, namely, q̃ = D_2g. Therefore, (Φ_2 q̃)^(2) = uv D_2(ρ v )+ vD_2(ρ uv ) + u D_2(ρ v^2 +p). Each term to the right of this equation reduces to zero for the given initial conditions. Since the right-hand-side of scheme <ref> is just a derivative of this, the scheme should sustain the initial condition.

In order to evaluate the one-dimensional shear flow numerically, we use Equations <ref> as initial conditions an a 64× 64 zone grid. We use an amplitude of the shear flow U=0.1, sound speed  c_s=1 on a domain of L_x,L_y=1. We run the simulation until t∼60, which corresponds to 60 sound crossing times. The results are shown in Figure <ref>, and it is evident that the scheme is maintaining the analytic solution. 







In this section, we investigate how the RBV2 algorithm maintains the one-dimensional shear flow presented in the previous section, with  a vertical hydrostatic equilibrium density gradient imposed. The equilibrium is set so that the pressure gradient balances an external gravity force, i.e., -dp/dy-ρ g=0. 

In the previous section the quantity q̃ is identically zero for the initial conditions. However, in the presence of a pressure gradient, q̃∼ D_2 p and therefore the RHS of Equation (<ref>) is not trivially zero. The analytic solution for the test case is,


    u(x,y)   =Usin(my) 
     
        v(x,y)   =0 
    
        p(x,y)   =c_s^2ρ(x,y) 
    ρ(x,y)   =ρ_0 exp(-gy/c_s^2)  .


We run the simulation using these as initial conditions, on a grid consisting of 64×64 zones. The simulation has dimensions L_x=L_y=1, c_s=1, ρ_0=1, and  g=1. The boundary conditions are ρ(x,y=0)=ρ_0, and ρ(x,y=1)=ρ_0/e, and are strictly periodic in the x direction. We show the results of the integration in Figure <ref>. 






§ SKEWED SHEAR FLOW


In this section, we investigate how the RBV2 scheme maintains a skewed shear flow, and present convergence tests to quantify the expected numerical dispersion in the algorithm.  We consider a flow velocity (u',v') in a  coordinate system (x',y'), rotated by some angle, θ, with respect to the unprimed coordinates.  The shear flow in the rotated system is  u'=A{e^iky'}. Therefore, in the unrotated (or computational grid) coordinates,


    u(x,y)   =Acos(θ){ e^ik(ycosθ-xsinθ)}
    
        v(x,y)   =Asin(θ) {e^ik(ycosθ-xsinθ)}
    ρ(x,y)   =ρ
    
        p(x,y)   =c_s^2ρ  .

We evolve these initial conditions on a 64× 64 grid, with L_x=L_y=√(2)π, and run the simulation for  60 sound crossing times, to be consistent with the one-dimensional shear flow analysis. We skew the shear flow by an angle θ=π/4. The results are shown in Figure <ref>. 



As shown in the residuals in Figure <ref>, there is  no dispersion in the skewed shear flow after 60 sound crossing times to machine precision. Assuming that the discrete vorticity equation has no dissipative terms,  if central difference schemes are used, then the numerical discretization does not produce
any dissipation of vorticity. In other words, only even-order derivatives are present when a truncation error analysis of the difference operators are made, or, one has a purely imaginary modified wavenumber.  However, the difference
operators do have dispersion error because their truncation error has odd-order derivatives and the
modified wavenumber for the first derivative deviates from being linear. This means that when compared to the analytic solution, there should be dispersion, and that the higher wavenumbers experience more rapid dispersion. 

In order to quantify the magnitude of this  numerical dispersion, we present a simulation of the same test with a higher wavenumber, which also may be thought of as a lower resolution simulation. In Figure <ref>, we increase the wavenumber by a factor of 4. We run this simulation also to 60 sound crossing times, and observe that as the dispersion manifests, the solution begins to go unstable, and the residuals grow to a factor of unity.




§ UPDATED SCHEME FOR EXTERNAL BODY FORCES


Many astrophysical applications require the advent of external body forces in numerical calculations. For example, a simulation of a protostellar disk requires the incorporation of both a coriolis deflection and an external gravitational field. However,  balancing source terms with non-zero flux gradients is known to be a difficult issue numerically <cit.>, even in the simplest case where the source terms cancel the flux gradients, such as a vertical hydrostatic equilibrium. We conducted numerical experiments that demonstrated that the RBV2 algorithm, when applied to a naive (cell-centered) source term discretization of an additional, constant gravitational force balanced by a pressure gradient, was unstable. Since the scheme was not developed to incorporate external body forces, it is possible to modify it to handle these and still retain minimal vorticity dissipation.

If a generic external acceleration a_ext=(0,a_x,a_y), was added to q̃, e.g.

    q̃=D_0μ w+D_1 f+D_2 g+μρ  a_ext  ,

the modified scheme  correctly incorporates the external force.  The external acceleration may be constant or spatially varying, and could possibly represent an external gravitational field. Furthermore, a_ext and ρ are defined on grid points, (j,k), so we include the μ operator so that each term in the residual is evaluated at the midpoints (j+1/2,k+1/2).

This new formula for the residual yields a straightforward inclusion of external forces in the scheme that may still be written as Λq̃ = 0. The proposed update requires that the difference and average operators commute.



We investigate how the updated scheme  maintains the one-dimensional shear flow presented in the previous section, in the presence of  a vertical hydrostatic equilibrium density gradient imposed. The equilibrium is set so that the pressure gradient balances an external gravity force (a_y=-g), i.e., -dp/dy-ρ g=0. The analytic solution for the test case is,


    u(x,y)   =Usin(my) 
     
        v(x,y)   =0 
    
        p(x,y)   =c_s^2ρ(x,y) 
    ρ(x,y)   =ρ_0 exp(-gy/c_s^2)  .


We run the simulation using these as initial conditions for ∼ 50 sound crossing times, on a grid consisting of 64×64 zones. The simulation has dimensions L_x=L_y=1, c_s=1, ρ_0=1, and  g=1. The boundary conditions are ρ(x,y=0)=ρ_0, and ρ(x,y=1)=ρ_0/e, and are strictly periodic in the x direction. We show the results of the integration in Figure <ref>. It is evident that the updated scheme maintains the analytic solution in the presence of the nonzero dissipation. 



§ NUMERICAL ASSESSMENT OF VORTICITY PRESERVATION FOR THE FULL EULER EQUATIONS




In this section, we numerically quantify the extent to which the RBV2 scheme preserves vorticity for the full Euler equations. We do this by examining the dissipation in the vorticity evolution equation for a spectrum of vortical modes up to the Nyquist frequency. For the full Euler equations, given by Equation <ref>, the scheme <ref> is vorticity preserving if (given by Definition 2 of <cit.>)


    D_0μ_0Ω̃+D_1[D_1μ(ρ uv)+D_2μ(ρ v^2)]-D_2[D_1μ(ρ u^2)+D_2μ(ρ uv)]=0  ,

where now Ω̃=D_1(ρ v)-D_2(ρ u). Therefore the vorticity equation is given by Equation <ref>, and the RHS of which  is RHS =h_1D_1curl_1+h_2D_2curl_2. Henceforth, we write

    q̃^(1)   =[D_0μ(ρ)+D_1(ρ u)+D_2( ρ v)]
    q̃^(2)   =[D_0μ(ρ u)+D_1(ρ u^2+p)+D_2 (ρ u v)]
    q̃^(3)   =[D_0μ(ρ v)+D_1(ρ uv)+D_2 (ρ v^2 +p)]  .


In this formalism, curl_1 = D_1(Φ_1 q̃)^(3)-D_2(Φ_1 q̃)^(2).Therefore (Φ_1 q̃)^(3)=uvq̃^(1)+vq̃^(2)+uq̃^(3) and (Φ_1 q̃)^(2)=(u^2+c_s^2)q̃^(1)+2uq̃^(2). Similarly, curl_2 = D_1(Φ_2 q̃)^(3)-D_2(Φ_2 q̃)^(2). To construct this,  (Φ_2 q̃)^(3)=(v^2+c_s^2)q̃^(1)+2vq̃^(3) and 
(Φ_2 q̃)^(2)=uvq̃^(1)+vq̃^(2)+uq̃^(3). Therefore


    curl_1=D_1[uvq̃^(1)+vq̃^(2)+uq̃^(3)]
        -D_2[(u^2+c_s^2)q̃^(1)+2uq̃^(2)]  ,


and


    curl_2=D_1[(v^2+c_s^2)q̃^(1)+2vq̃^(3)] 
        -D_2[uvq̃^(1)+vq̃^(2)+uq̃^(3)]  .


Then, for a uniform grid where h_1=h_2=h the RHS simplifies to 


    RHS=h[D_1D_1-D_2D_2][uvq̃^(1)+vq̃^(2)+uq̃^(3)]
    +hD_1D_2[(v^2-u^2)q̃^(1)-2uq̃^(2)+2vq̃^(3)]  ,

with μ_0=μ^2, μ=μ_1μ_2, D_0=1/Δ tΔ̅, D_1=1/δ xδ_1μ_2 and D_2=1/δ yδ_2μ_1, as defined in Section 2. If the dissipation in the RBV2 scheme exactly preservers vorticity, Equation <ref> should yield identically zero on the grid. Any deviations would signify that the scheme was dissipating vorticity. 

We calculate the residuals of the RHS from Equation <ref>. The results of this are shown in Figure <ref>. Using a 32×32 zone grid, we seed vortical perturbations of the form



    u(x,y)   =Ak_x/√(k_x^2+k_y^2){exp(i2π/L(k_xx-k_yy))}
     
        v(x,y)   =Ak_y/√(k_x^2+k_y^2){exp(i2π/L(k_xx-k_yy))}
    
        p(x,y)   =c_s^2ρ(x,y) 
    ρ(x,y)   =ρ_0    ,

for a set of all integers (k_x,k_y)∈ℤ×ℤ from 1 to the nyquist wavenumber N/2=16. We evaluate the RHS of the discrete vorticity equation, and show the  root mean square of this quantity normalized by the initial vorticity, <RHS>/<Ω̃>^2 in Figure <ref>. We show the linear and log scale of this assessment, and in the linear scaling include the case where k_x=k_y=0, which correspond to the unskewed shear layers that we have analytically shown to be zero. It appears that RBV2 perfectly evolves vorticity for vortical modes of equal wavenumbers, e.g. k_x=k_y. The scheme dissipates vorticity when the wavenumbers are antisymmetric, e.g.k_x k_y , and the dissipation increases for higher wavenumbers. 



§ DYNAMICAL INTERACTION OF VORTICES IN A PROTOPLANETARY DISK



There now exists mounting theoretical and observational evidence that vortices play  an important role in the local and global dynamics of protoplanetary disks. Therefore, developing computational algorithms that can accurately evolve flows replete with vortices and turbulence is critical to any theoretical understanding of them.  Vortices have been invoked as the sites for  planet formation due to the rapid settling of dust grains into their central regions <cit.>. This process enhances the local surface density, and may trigger giant planet formation via gravitational instability <cit.>.<cit.> reported ALMA observations of the disk around the star  Oph IRS 48 that traced millimeter sized grains, and were consistent with a massive, anticyclonic dust trapping vortex. <cit.> presented similar ALMA observations as evidence for a massive (≥0.3M_Jup) vortex in the HD 135344B disk. Moreover, they argued that  the  vortex may be massive enough to launch spiral density waves, and could therefore be responsible for  the observed, global, spiral structure.  <cit.> suggested that vortices in hot accretion disks may explain the variability in their X-ray temporal power spectra. Similar temporal variability is observed in young stellar objects <cit.>, which may also be caused by vortices. 

From a theoretical standpoint, vortices and hydrodynamical turbulence in protostellar disks may be important drivers of  angular momentum transport. In the ideal MHD limit, the magnetorotational instability <cit.> produces both the  magnitude and outward directionality  of the angular momentum transport necessary for the viscous evolution of  gas in an accretion disk. However, protoplanetary disks are so weakly ionized that magnetic fields may not affect the motion of the neutral medium, so there remains interest in searching for a purely hydrodynamical mechanism for the outward transport of angular momentum.  <cit.> discovered that vortices that form naturally above the midplane from disk perturbations  are robust for hundreds of orbit, and could produce significant  angular momentum transport. Self-gravitating, cooling instabilites <cit.>, baroclinic instabilities <cit.> and zombie vortex instabilities <cit.> may also transport angular momentum (for a recent review, see <cit.>). 
 
Here, we present a simulation of the dynamical interaction of seven vortices in a protoplanetary disk. This provides an example of how the RBV2 algorithm, updated to handle external body forces, will be useful for studying turbulent astrophysical flows. We perform a simulation in the Cartesian shearing box formalism, which is employed to examine local phenomenon at a scale which captures the shear, yet omits the global curvilinear geometry, described in detail in <cit.>. For a Keplerian disk, the hydrodynamical equations of motions expanded around a fiducial disk radius are 



    ∂ v/∂ t +  v·∇⃗ v=-∇⃗ P/ρ-2Ω× v +3|Ω|^2 x x̂  ,

where v, P, ρ and Ω are the local velocity, pressure, density and rotational velocity. In addition to the descripton of the implentation in Section 4.1 in <cit.>, the operator q̃ has been updated at each grid cell to include the appropriate coriolis term, -2Ω× v, and   tidal expansion  of the effective centrifugal and gravitational potential from the central object term, 3|Ω|^2 x x̂. The seeded vortices have the functional form proposed in <cit.>, and used in Section 4.2 and 4.3 of <cit.>. The simulation  was run on 512×512 zones, at a fiducial radius of 1AU, with a sound speed c_s=0.05 and radial and azimuthal extent of .04AU. 

Figure <ref> shows the dynamical interaction of  seven anticyclonic vortices in a protoplanetary disk as rendered with the updated RBV2 algorithm. Only anticyclones are expected to be long-lived in disks, because they have the same sense of rotation as that of the local shear, which is opposite that of the global rotation. The same shear which enforces the rotation will rip apart  cyclonic vortices <cit.>. In the absence of the background shear and the additional forces present in the rotating frame, vortices  that have the same spin orientation and are close enough together will eventually merge <cit.>. As demonstrated in the second and third panels of Figure <ref>, the five central vortices dynamically interact very rapidly. The Keplerian shear  dominates the dynamics of the two furthest vortices, driving the upper left vortex upwards and lower right vortex downwards. In the fourth panel, the five central vortices have coalesced into two new vortices with substantially more sub-structure. Due to the periodic boundary condition in the azimuthal direction, the two furthest vortices have also begun merging at this stage. The three long-lived vortices continue to shear past each other for ∼ 10 orbits, dynamically interacting, exchanging vorticity and shedding Rossby waves that dramatically effect the sub-structure in the surrounding flow. After ∼ 18 sound crossing times, the three vortices coalesce into a single, much larger whirlpool.






§ CONCLUSION


<cit.> presented a novel approach to the accurate evolution of vorticity by insisting on vorticity preservation in idealized hydrodynamical cases on a multi-dimensional grid. In this paper, we investigate the vorticity-preserving properties of the resultant RBV2 algorithm. We evaluate how the scheme handles several idealized hydrodynamical tests of the full Euler equations. We propose an adjustment to the dissipation that retains the vorticity-preserving qualities, while accurately incorporating external body forces, and demonstrate that the adjusted scheme sustains a hydrostatic equilibrium. 
We present a novel numerical assessment of vorticity dissipation for discrete wave-number, vortical modes. We find that for the full Euler equations, RBV2 perfectly preserves vorticty for modes with symmetric wavenumbers, and that the error increases with asymmetry. We then perfrom shearing sheet simulations of vortex interactions in a protoplanetary disk, to demonstrate the utility of the updated algorithm for modeling astrophysical fluid flows. 

The RBV2 algorithm has proven to exhibit remarkably efficient and accurate results for test problems involving turbulence and vortices. Given its minimal dissipation of vorticity, we plan use  the  algorithm to simulate astrophysical environments  where turbulence is expected to dominate the dynamics. One very promising example of this is a protoplanetary disk, where turbulence and vortices may dictate the accretion process <cit.>. With the advent of the proposed update to include body forces, we believe that the  scheme will prove extremely useful for elucidating the nature of and interaction between vortices in Jupiter's polar region, where the small Rossby number <cit.> implies that the Coriolis deflection dominates the dynamics.


This material is based upon work supported by the National
Aeronautics and Space Administration through the NASA
Astrobiology Institute under Cooperative Agreement Notice
NNH13ZDA017C issued through the Science Mission Directorate.
We acknowledge support from the NASA Astrobiology
Institute through a cooperative agreement between NASA
Ames Research Center and Yale University.

We thank Professor Alain Lerat for extensive correspondence regarding the contents of this paper. We thank Professor Gregory Laughlin for extensive support and insightful conversations. 

https://github.com/DSeligman/MAELSTROM2D




































































