Multidimensional, compressible viscous flow on a moving Voronoi mesh
Abstract
Numerous formulations of finite volume schemes for the Euler and NavierStokes equations exist, but in the majority of cases they have been developed for structured and stationary meshes. In many applications, more flexible mesh geometries that can dynamically adjust to the problem at hand and move with the flow in a (quasi) Lagrangian fashion would, however, be highly desirable, as this can allow a significant reduction of advection errors and an accurate realization of curved and moving boundary conditions. Here we describe a novel formulation of viscous continuum hydrodynamics that solves the equations of motion on a Voronoi mesh created by a set of meshgenerating points. The points can move in an arbitrary manner, but the most natural motion is that given by the fluid velocity itself, such that the mesh dynamically adjusts to the flow. Owing to the mathematical properties of the Voronoi tessellation, pathological meshtwisting effects are avoided. Our implementation considers the full NavierStokes equations and has been realized in the AREPO code both in 2D and 3D. We propose a new approach to compute accurate viscous fluxes for a dynamic Voronoi mesh, and use this to formulate a finite volume solver of the NavierStokes equations. Through a number of test problems, including circular Couette flow and flow past a cylindrical obstacle, we show that our new scheme combines good accuracy with geometric flexibility, and hence promises to be competitive with other highly refined Eulerian methods. This will in particular allow astrophysical applications of the AREPO code where physical viscosity is important, such as in the hot plasma in galaxy clusters, or for viscous accretion disk models.
keywords:
hydrodynamics – methods: numerical.1 Introduction
The last two decades have seen remarkable advances in the numerical solution of the compressible NavierStokes (NS) equations, which lies at the heart of computational fluid dynamics (CFD) and computational aeroacoustics, but also as numerous applications in astrophysics. In particular, important progress has been made in approaches based on the finite volume method (FVM), both using structured as well as unstructured grids (see Mavriplis, 1997, for a review). Other popular techniques include finite element methods (FEM), discontinuous Galerkin schemes, and even meshfree approaches such as smoothed particle hydrodynamics (Sijacki & Springel, 2006).
When unstructured grids have been employed, they were most most often in the form of triangular grids in two dimensions, or tetrahedral grids in three dimensions. Indeed, finitevolume implementations of the twodimensional NS equations on triangular meshes date back to work by Mavriplis & Jameson (1990), Frink (1994) and Coirier & Powell (1996). Much recent work has also focused on developing optimum meshgenerating algorithms that require minimal human input and yield efficient representations of geometrically complex simulation domains. However, little work has been done on dynamically evolving meshes, such as those we shall consider here.
Because unstructured meshes have been demonstrated to be accurate and efficient for both steadystate and transient compressible inviscid flows (Barth, 1992; Venkatakrishnan, 1996), they are now used regularly in engineering applications. Moreover, the geometric flexibility of unstructured grids allows the use of simple coordinate systems (in the laboratory frame) without the need to work with complex coordinate transformations to describe curved surfaces (e.g. see Toro, 2009). Indeed, hard boundaries can be tailored by carefully positioning a few cell faces or mesh generating points along the surface, and creating the triangulation through Delaunay tessellation. As a result, most NS applications on unstructured meshes for industrial design make use of triangular grids, typically based on the finite element method, although finite volume schemes have also been considered. Detailed reviews and stability analysis of explicit FVM for the NS equations on Cartesian and Delaunay meshes can be found, e.g, in the doctoral theses of Coirier (1994) and Munikrishna (2009).
In this work, we present a numerical scheme that solves the NS equations on a general unstructured moving mesh that is constructed as the Voronoi tessellation of a distributed set of points that move with the local velocity field. Despite being, in the general sense, an “unstructured” mesh, the Voronoi diagram has a mathematically welldefined structure that makes the resulting schemes comparatively simple and robust (e.g. Mishev, 1998). In fact, this type of mesh is commonly adapted for the construction of finite volume methods for elliptic problems and has been in use in numerical studies of solid state physics (Sukumar et al., 1998; Sukumar, 2009) such as simulations of fractures and cracks (Sukumar & Bolander, 2009), as well as numerical simulations of oil reservoirs. Some studies (Christov, 2009) have also examined how reconstructions designed for unstructured triangulations can be extended to static Voronoi meshes.
However, Voronoi meshes have infrequently been applied to hyperbolic conservation laws such as the Euler equations, let alone moving Voronoi meshes. To our knowledge, the earliest attempts to use dynamically adaptive Voronoi tessellations for the NS equations data back to Borgers & Peskin (1987), although for very simplified, incompressible, twodimensional problems. Around the same time, Dukowicz et al. (1989) developed the General Topology Godunov Method. This method – based on a mesh that is not quite a Voronoi tessellation, but similar in spirit – was introduced as an alternative to the Lagrangian particle methods (see, for example Brackbill & Monaghan, 1988) which gained increasing popularity in computational plasma physics and astrophysics in the following years.
Recently, a complete threedimensional implementation of the Euler equations on a moving Voronoi mesh has been described and implemented in the AREPO code by Springel (2010) (see also Duffell & MacFadyen, 2011). The work we present here is an extension of the AREPO scheme to the NS equations, which we have realized in this code as an optional module. AREPO can be classified as an arbitrary Lagrangian/Eulerian (ALE; Hirt et al., 1974) code, in the sense that the mesh can be moved with the velocity of the flow so that quasiLagrangian behavior results and the mass flux between cells is minimized (although it is not strictly zero, in general). On the other hand, the mesh may also be kept stationary if desired, effectively yielding an Eulerian formulation. We note that because the meshgenerating points may also be arranged on a regular lattice and arbitrarily refined with time, the AREPO code naturally includes ordinary Eulerian techniques on a Cartesian grid and adaptive mesh refinement (AMR) algorithms as special cases.
Besides the work of Duffell & MacFadyen (2011), the new VoronoiALE method of Norris et al. (2010), which includes viscous terms, is the approach most closely related to that presented here, although it is restricted to the incompressible NS equations. Also, Ata et al. (2009) have applied a Voronoibased finite volume scheme to the twodimensional inviscid shallow water equations, in terms of an algorithm they referred to as the ‘natural volume’ method.
Although primarily designed for astrophysical fluid dynamics where selfgravity is an important ingredient (see for example Vogelsberger et al., 2011), the moving Voronoi mesh approach of AREPO offers a number of features than can be advantageous for more general problems in fluid dynamics. First, the moving mesh geometry is adaptive in a continuous manner and can naturally respond to the local flow, increasing the resolution automatically and smoothly in regions where the flow converges. (In contrast, AMR codes refine the grid discontinuously in time, which can introduce errors that are potentially difficult to assess.) Importantly, this Lagrangian character of the dynamics yields reduced advection errors and a very low numerical diffusivity of the scheme. Second, the moving mesh formulation retains the Galileaninvariance of the fluid dynamics at the discretized level of the equations (Springel, 2010). In other words, the truncation error of the scheme does not depend on the bulk velocity of the system, unlike for traditional Eulerian and AMR codes, and the quality of the solution does not degrade when highspeed flows are present. While conventional fixedmesh Eulerian codes may, in principle, be able to suppress additional errors from large bulk velocities by using a sufficiently fine mesh (see Robertson et al., 2010, for a study of Galilean invariance in grid codes), this strategy can become computationally prohibitive, and it also depends on the magnitude of the bulk velocity involved. It is therefore desirable to construct efficient methods that yield manifestly Galileaninvariant solutions (modulo floating point roundoff errors). Third, the moving mesh approach allows much larger timesteps in the case of rapidly moving flows, because it can avoid the stability constraint (where is the cell size and the bulk velocity) that augments the Courant condition in the Eulerian case.
From an astrophysical standpoint, compressible viscous flow remains a viable approximation to more complex or computationally expensive momentum transport mechanisms such as magnetohydrodynamic turbulence or anisotropic plasma viscosity. Global simulations of cold accretion disks around protostellar objects (e.g. see de ValBorro et al., 2006) still include shear viscosity coefficients in the form of a ShakuraSunyaev eddy viscosity coefficient (Shakura & Sunyaev, 1973).
An even clearer case for the need of a viscous treatement of astrophysical gasdynamics is given by the interacluster medium of hot galaxy clusters. Here the SpitzerBraginskii viscosity (Braginskii, 1965) becomes quite significant, certainly in the unmagnetized case, which has been studied both using grid (Ruszkowski et al., 2004) and SPH (Sijacki & Springel, 2006) codes. In this regime, the commonly adopted assumption of inviscid behaviour with an effectively infinite Reynolds number is in principle incorrect and should in future simulation work be replaced with a full accounting of the correct physical viscosity.
Additionally, physical viscosity can be implemented on turbulent cascades with resolved inertial range (see Bauer & Springel, 2011, for an application of our viscosity approach) in order to prescribe a wellspecified Reynolds number and a physically correct shape for the dissipation range, unaffected by the details of the numerical viscosity of the hydro scheme, which would otherwise induce the dissipation of turbulence on the grid scale. This can in particular inform the ongoing debate whether artificial viscosity effects in SPH can affect the turbulent cascade above the formal resolution length (Bauer & Springel, 2011; Price, 2012).
This paper is organized as follows. In Section 2, we briefly review the basic NS equations we want to solve, and the role and meaning of the different viscosity coefficients. In Section 3, we then introduce in detail our discretization and time integration schemes, emphasizing a description of the calculation of suitable velocity gradient estimates at face centers, and of secondorder derivatives of the velocity field. We then move on to discuss the performance of our new approach for a number of test problems in Section 3. Finally, we summarize our results and present our conclusions in Section 4.
2 The NavierStokes Equations
The compact form of the Euler equations, when written in terms of the vector of conserved quantities (Toro, 2009) is
(1) 
with
(2) 
and where
(3) 
is the massmomentumenergy flux density tensor (). The operator in Eq. (1) is a tensor divergence, i.e. in tensor notation we have . The momentum components in the conservative form of Equation (1) represent a transfer of momentum, owing merely to the mechanical transport of different particles of fluid from place to place and to the pressure forces acting on the fluid (e.g. Landau & Lifshitz, 1959). In Eq. (1) we have made the advective character of the fluxes explicit by denoting them .
The internal friction present in any real fluid causes an irreversible transfer of momentum from points where the velocity is large to those where it is small. The momentum flux density tensor is thus altered from its ideal from in Eq. (3), where it only contains an inertial and an isotropic component (described by a symmetric stress tensor due to the local pressure ), to a modified expression that accounts for an irreversible viscous transfer of momentum
(4) 
where is the total stress tensor and is called the viscous stress tensor. The latter includes the effects of isotropic compression and expansion forces (“bulk viscosity”) as well as shearing forces (“shear viscosity”).
Similarly, the energy component of Eq. (3) is affected by the inclusion of the viscous stress tensor. Because of the dissipative nature of viscosity, a conservative formulation of the NS equations must include a contribution of to the energy budget, i.e. the work per unit area per unit time,
(5) 
needs to explicitly account for the work done by viscous forces.
A general parametrization of the viscous stress tensor is given by
(6) 
Often, the viscous stress tensor is decomposed into a traceless part and a diagonal part, such that the first corresponds to constantvolume shear deformations (often called the rateofdeformation tensor) and the second to isotropic expansions/contractions. Accordingly, in Eq. (6) is commonly referred to as the shear viscosity and as the bulk viscosity. The degree of resistance to uniform contractions/expansions is intrinsic to the molecular/chemical properties of the fluid in question, and can be understood through kinetic theory. In this picture, bulk viscosity arises because kinetic energy of molecules is transferred to internal degrees of freedom. Ideal monoatomic gases (modeled as hard spheres interacting only through elastic collisions) have no internal degrees of freedom, and are thus expected to have vanishing bulk viscosity. At one time Stokes suggested that this might in general be true (the socalled Stokes’ hypothesis of ) but later wrote that he never put much faith in this relationship (Graebel, 2007). Indeed, when deviations from the ideal gas equation of state are included in a hardsphere, ChapmanEnskog approach to kinetic theory, a nonzero value for the bulk viscosity is obtained. In an extension of the hard sphere fluid model, the LonguetHigginsPople relation results (March, 2002), motivating the hypothesis that both viscosities are always related in a linear fashion (but see Meier et al., 2005). In general, we consider and as essentially arbitrary input properties to our simulations, which may also depend on local physical parameters such as temperature or density. Although the effects of physical bulk viscosity are not harder to implement numerically than those of shear viscosity, the physical origin of bulk viscosity is often less clear. Also, we note that many numerical solvers for viscous flow focus on the incompressible regime (), where the existence of a physical bulk viscosity is in any case not of importance. However, for compressible flow, the value of may still become important in certain situations.
When the effects of viscosity are included, the formerly homogeneous differential equations of the Euler form (Eq. 1) become
(7) 
where is a viscous source term given by
(8) 
The solution of the Euler equations with source terms is often handled by operatorsplitting methods (e.g. Toro, 2009; LeVeque, 2002). That is, the numerical scheme alternates between an advective step that solves the homogeneous part, and a sourceterm step. Thus, the solution of Eq. (7) is split into a two stage problem:
(9)  
(10) 
Typically, the source terms are more easily written in the primitive variable formulation of the Euler equations. A common choice of the primitivevariable vector is , which we also adopt here. For sources corresponding to the NS viscous terms (Eq. 8), only the component of is affected, thus simplifying the solution method of the sourceterm step. The threedimensional Euler equations can be written in the primitive variable form as (Toro, 2009)
(11) 
For this choice of variables, the coefficient matrices are given by (Toro, 2009)
(12)  
(13)  
(14) 
which is exactly equivalent to the familiar equations
(15a)  
(15b)  
(15c) 
In this formulation, the viscous terms of the NS equations, which affect only the velocity, are (e.g. Landau & Lifshitz, 1959)
(16) 
An alternative to expressing the viscosity effects as source terms is to absorb them directly into the flux divergence,
(17) 
which highlights the still conservative character of the NS equations. Here diffusive fluxes, defined by
(18) 
are responsible for the effects of viscosity. An implementation of the diffusive fluxes in this conservationlaw form is clearly the preferred choice for FVM schemes, which are specifically designed for solving the integral form of these conservation laws. In fact, in this case they exactly conserve all the involved quantities to machine precision. We will therefore focus on this method in our study. The central aspect will be the numerical scheme used for estimating the velocity gradients at the cell interfaces, and hence the discretization of the diffusive fluxes. In the next section, we describe our approach for this in detail.
3 A Finite Volume Scheme with Viscous Fluxes on a Voronoi Mesh
3.1 Basic MUSCLHancock Finite Volume Scheme: Overview
Finite volume methods enforce the integral form of the conservation laws on discrete meshes. This approach is manifestly conservative, since fluxes of quantities that leave a cell simply enter the neighboring cell. The NS equations in finitevolume form are
with  (19) 
where, in general, the intercell fluxes contain both advective and diffusive contributions,
(20) 
The scheme used by AREPO is the finite volume MUSCLHancock approach, consisting of a MUSCL (Monotone Upstreamcentered Schemes for Conservation Laws) linear reconstruction stage, and a Hancock twostage time integration
(21) 
where the numerical fluxes represent appropriately timeaveraged approximations to the true flux across the interface shared by cells and . The time label in Eq. (21) indicates that an intermediatestage (a half timestep evolution) has been performed to obtain the numerical estimate of , meaning that the timestepping in Eq. (21) uses timecentered fluxes, giving it secondorder accuracy. The Hancock part of the scheme is a twostep approach (the familiar predictorcorrector algorithm) in which the correction halfstep is obtained from the solution of the 1D Riemann problem across each face of the control volume. The general finite volume MUSCLHancock scheme has hence the following three steps (Toro, 2009):
(I) Gradient Estimation, Linear Data Reconstruction and Boundary Value Extrapolation
Once a local gradient estimate for the conserved quantities of cell is available, linear data reconstruction takes the form
(22) 
where we denote by the estimated vector of conserved variables at the centroid of the interface, obtained by linearly extrapolating the cellcentered values of the th cell (on the “left” side) from , the cell’s center position, to . Similarly, corresponds to the estimates of the facecentroid values obtained by linear extrapolation of the cellcentered values of the th cell (the “right” side), whose center position is . In both cases, is the position vector of the face centroid between the cells. The Jacobian is explicitly labeled with superscript to point out that it corresponds to the estimate of spatial derivatives at the beginning of the timestep.
(II) Evolution of Boundary Extrapolated Values
This is, strictly speaking, the “predictor” half timestep. The conserved variables are evolved for with flux estimates obtained from the values at the beginning of the timestep:
(23) 
(III) Solution of 1D Riemann Problems and Computation of Godunov Fluxes
This corresponds to the “corrector” half timestep in the twostage Hancock approach. Once the values to the right and left of the interface at time are known, the discontinuity is treated as a onedimensional Riemann problem. An exact or approximate Riemann solver is used to return values of , and at the interface, at a time corresponding to . From these values, the advective fluxes can be directly computed (Eq. 3). These are timecentered fluxes used to update the system from the beginning of the timestep to its end,
(24) 
Figures 1 and 2 illustrate the mesh geometry and the basic steps of this inviscid numerical scheme implemented in AREPO. One additional point we have not explicitly discussed here for simplicity is the treatment of the mesh motion, as indicated in Fig. 2. This is incorporated into the scheme by evaluating all fluxes in the rest frame of the corresponding face, as described by Springel (2010). This requires appropriate boosts of the fluid states and the fluxes from the lab frame to the rest frame of each face, and back. For a Voronoi mesh, the face velocities are fully specified by the velocities of all the mesh generating points. The latter can be chosen freely in principle, but if they are set equal to the fluid velocities of the corresponding cells, a Lagrangian behavior and a manifestly Galileaninvariant discretization scheme is obtained in which the truncation error does not depend on the bulk velocity of the system.
3.2 A MUSCLHancock FiniteVolume Scheme with Viscous Terms
A cellcentered, finitevolume solution of the NS equation can be written as
(25) 
where we have retained the distinction between advective and viscous fluxes. As in the case of the Euler equations, the numerical method essentially consists of the problem of finding accurate timecentered numerical fluxes across each of the interfaces of a cell. How to do this in detail for the diffusive part of the fluxes has been the focus of numerous efficiency and stability analyses (see Puigt et al., 2010, for a detailed description).
Eq. (25) uses timecentered fluxes, obtained here with the twostep Hancock technique, as described above. Thus, for estimating both and a half timestep predictor stage is required. In the MUSCLHancock approach for inviscid flow, this step is carried out by linear reconstruction from each cell center to the interface, followed by solving a onedimensional Riemann problem at the interface where the extrapolations meet. The traditional formulation of the Riemann problem and its solution are exclusive to hyperbolic differential equations and thus do not provide exact solutions for the NS equations. Since a general solution for the viscous Riemann problem does not exist, we will treat the viscous fluxes in Eq. (25) as a correction to the solution of an otherwise inviscid flow.
Our NS version of the MUSCLHancock scheme consists of the following three different stages (in addition to those described in Section. 3.1):

Correct the MUSCL linear extrapolation of primitive variables by applying a viscous kick.

Extrapolate the cellcentered gradients linearly and evolve them for half a timestep.

Average the extrapolated velocity gradients at the interface and use them to estimate viscous fluxes.
To extrapolate the gradients from their cellcentered values to the interfaces, information about the higherorder derivatives of the primitive variables is needed. If gradients are assumed to vary linearly in space, an estimator for the Hessian matrix for each of the five primitive variables is sufficient. Evidently, enough information is contained in the cellcentered quantities to estimate both the local gradient and the Hessian corresponding to a given scalar quantity . However, estimating both of these simultaneously is significantly more difficult than estimating them one after the other. Therefore, we will effectively treat and as two independent fields that vary linearly in space, and this variation needs to be estimated from the mesh data through a suitably discretized differential operator.
As a simpler alternative to the gradient reconstruction approach, we briefly describe how one can use the gradients already available from the linear reconstruction step. In this approximation, a given quantity varies only linearly within the control volume, such that consistently evaluated gradients are piecewise constant. This means that each interface represents a discontinuity in the gradient field . Naively, one may think that the arithmetic average of both gradients that meet at a face is a good estimate for the gradient at the interface itself. However, on second thought, one realizes that both cells do no necessarily have the same weight if cells of different volume meet. Furthermore, the unweighted average of the two cellcentered values really represents the value at the midpoint of the two meshgenerating points, which, for a Voronoi mesh, can be substantially offset from the midpoint of the face. We therefore adopt the approach of Loh (2007), which consists in choosing one of the two gradients that meet at the interface, based on prior knowledge of the direction of the flow across the interface. Thus the threestage scheme introduced above could be alternatively replaced by the simpler method:

At the cell interface where two different gradients meet, choose the upwind gradient.
In either method, once we have an estimate of both viscous and advective fluxes, the timestep evolution of the conserved quantities is carried out as in Eq. (25). However, the approach (AC) is preferable to the Loh (2007) scheme because it uses timecentered estimates for both and , hence preserving the order of accuracy of the original inviscid scheme. We therefore now provide a more detailed description of the individual steps in this threestage approach.
(A) Viscosity Kicks
Although Eq. (25) is written in an unsplit form, the predictor step is indeed operator split, evolving the advective and diffusive terms separately (e.g. Coirier & Powell, 1996). While our method for estimating the advective fluxes remains the MUSCLHancock scheme, the technique for estimating the diffusive fluxes is essentially contained in the estimation of the velocity gradients at each interface (see Coirier, 1994; Puigt et al., 2010, for a series of tests on different interface gradient estimates). Looking for better accuracy, we have chosen to couple these two otherwise independent procedures by correcting/biasing the linear extrapolation of the velocity field (stage in Section 3.1) with a viscous source term.
The benefit of carrying out a linear extrapolation to cell interfaces in primitive variables is the simplicity of the Galilean transformation needed to boost the quantities to the frame of a moving interface. Since the Galilean boost does not affect the mass and pressure of a given cell, only the local velocity field is transformed. In addition, adding force source terms to the equations of motion in primitive variable formulation is simpler, since these only couple to the momentum equations. Thus, a “viscous kick” can be applied to the velocity field in the half timestep evolution stage:
(26) 
In this way, the subsequent linear extrapolation of primitive variables will already include viscosity effects to first order in time.
While working with numerical fluxes across interfaces requires velocity gradients, the use of cellcentered source terms in Eq. (26) calls for second order derivatives of the velocity field. Thus, in addition to the cellcentered velocity gradients , and , the cellcentered Hessian matrices , and are now needed. As we will see below, these matrices will be of use in more than one occasion, justifying the computational cost incurred to calculate them.
(B) Linear Extrapolation of Gradients
The linear reconstruction implemented in our MUSCLHancock approach essentially assumes that the gradient of a scalar quantity does not vary significantly across the spatial scale of a cell. For smooth flows, the gradients of two neighboring cells and will not differ significantly. Furthermore, in the presence of strong discontinuities, gradients on each side will be slopelimited, and therefore will not differ by much either. Hence, a first guess for the gradient at the interface between two cells is just the average of the cellcentered estimates at each side of the face
(27) 
However, as we pointed out earlier, the gradient average above is actually representative of the midpoint between the two cell centers and , which in general does not lie close to the center of the face in a Voronoi mesh, and may in fact lie within a third cell. Unless gradients are assumed to vary within a cell, it will not be possible to assign the estimate to the center of the interface with any confidence.
Let us assume that the scalar field is infinitely differentiable and, consequently, so is its first derivative. Thus, we can Taylor expand both quantities to arbitrary order around a mesh generating point :
(28)  
(29) 
where is the Hessian matrix of the scalar quantity and is a tensor containing the thirdorder derivatives of (i.e. ). Truncating both Taylor expansions to firstorder in , we see that we can obtain linear reconstructions for both the physical quantities and their gradients provided that we have numerical estimates for both the gradients and the Hessians at each mesh generating point.
We emphasize that a Taylor expansion is not equivalent to a polynomial data reconstruction. Indeed, it is desirable that reconstruction schemes are manifestly conservative, in the sense that the average of the reconstruction over the cell should be identical to the value of at the geometric center of the cell. This property of reconstruction schemes is sometimes referred to as exactness, meaning that if a polynomial reconstruction is cellaveraged over the mesh, the reconstruction procedure recovers the same polynomial. This condition is trivially satisfied for a linear reconstruction of the form . However, higherorder reconstruction schemes require the use of zeromean polynomials, which, beyond firstorder, differ from the Taylor series (e.g. Colella & Woodward, 1984; Coirier & Powell, 1996).
The linear reconstruction of the scalar field and of the vector field , treated as if they were independent quantities, effectively constitutes a hybrid method between standard linear reconstruction and fully Kexact secondorder reconstruction, as illustrated in Figure 3. In this approximation, second derivatives are considered negligible for the spatial reconstruction of the primitive quantities, but they are still included for a more accurate estimate of the gradients near the cell interfaces. We also note, that in this way our numerical scheme reduces to that originally in AREPO (which is secondorderaccurate) when the viscous fluxes are disabled.
Once an estimate for the Hessian matrix is available (Section 3.3), a linear extrapolation of the gradients from the cell centers to the interfaces can be obtained from
(30) 
which is a better approximation than Eq. (27). However, the time evolution of the gradients during a single step could be equally important as their spatial variation over the length scale of a cell, hence we also need to evolve them for half a timestep to obtain a time integration scheme that is consistent with the secondorder accurate twostage MUSCLHancock approach. In the latter, to extrapolate and evolve a scalar quantity we consider
(31) 
where the time derivative of the quantity in the control volume of the th cell can be obtained from the primitive variable formulation of the Euler equations in tensor notation:
(32) 
Here sums over repeated indices are implied. Latin indices take the values or , while Greek indices take the values and are used to number the components of the primitive quantity vector ( for , respectively). As with our previous notation, the indices and are reserved for labeling the mesh generating points and their associated cells.
Eq. (32) is an advection equation for the primitive variables. Analogously, to “advect” the gradients of the primitive variables from the cell center to the interface, we can ignore the viscous terms and derive an equation of motion for the spatial derivatives by differentiating Eq. (32):
(33) 
where we can identify the Jacobian matrix of the primitive variables as , and the Hessian tensor () of the primitive variables as . Therefore, the time derivative of each component of the primitive variable Jacobian matrix is
(34) 
where we introduced the rank tensor . Since is a function of the primitive variables , the tensor can also be written as (see Appendix)
(35) 
and therefore its numerical estimate is given by the product of the exact derivatives (evaluated with values of the primitive variables at the center of the cell) and the (already available) numerical estimates for the gradients . The second term on the right hand side of Eq. (34) is the product of the known coefficients (evaluated at the center of the cell) and the numerical estimates of the Hessian tensor .
Finally, with a numerical estimate of at hand (see Section 3.3), the extrapolated and half timestep evolved gradients of the velocity are (in analogy to Eq. 31):
(36) 
with analogous expressions for and . In Eq. (36), the term is obtained from Eq. (34) with and .
In Fig. 4, we show a sketch of the different steps involved in obtaining timecentered diffusive fluxes. We point out that taking the Hessian matrices of the velocity field to be identically zero is equivalent to the alternative scheme . The third term to the right hand side of Eq. (36) is still different from zero even if (Eq. 34) since, in general, . By advecting the gradients according to Eq. (34) we gain additional accuracy at no additional computational expense because the terms are known exactly (see Appendix), given the values of the primitive variables and their respective gradients at the center of each cell.
(C) Viscous Flux Calculation
An accurate estimate of the viscous fluxes between two cells requires an accurate estimate of the velocity gradients at the interface. The gradient extrapolation method described above produces in general two different values of the velocity gradient that meet at the interface. This defines a general Riemann problem for the differential equation in Eq. (34) which is no longer a homogeneous hyperbolic differential equation. Therefore, attempting to solve this new Riemann problem for the spatial derivatives of the scalar quantities introduces a significant additional difficulty. For simplicity, we will assume that the differences between two gradient extrapolations meeting at an interface are small enough such that a simple arithmetic mean can be used. This assumption, of course, is valid only when the field of second derivatives is sufficiently smooth (see Section 3.3).
The time and area averaged flux across the face  that moves with speed is defined as
(37) 
The numerical or Godunov estimate of these fluxes is chosen so that the analytic expressions for and are evaluated with numerical estimates of , and at the centroid of the interface. The advective Godunov fluxes are
(38) 
where is the conserved variable vector at the centroid of the interface, as seen in the lab frame, obtained from the solution of a 1D Riemann problem across the  interface and along its normal. Multiplying by is equivalent to projecting the flux matrix (Eq. 3) along the normal of each face. The Godunov fluxes and are thus component vectors. The diffusive Godunov flux vector is obtained from the diffusive flux matrix
(39) 
where are the components of the viscous stress tensor , which depend on the local value of the velocity and the velocity gradients. These components are:
(40) 
Just like with the advective fluxes, the flux tensor (Eq. 39) must be projected onto the normal of each interface to obtain the 5component vector
(41) 
where is the primitive variable vector at the centroid of the interface, as seen in the lab frame (whose associated conserved variables are in Eq. 38). The spatial derivatives correspond to our extrapolateandaverage scheme for linearly varying gradients. As with , we are interested in estimates of at the centroid of the face. For both these quantities, only the velocity and its spatial derivatives are relevant when viscous fluxes are calculated.
3.3 Hessian Estimation
In analogy to the gradient calculation for Voronoi meshes discussed by Springel (2010), here we discuss the estimates of the cellcentered Hessian matrices for each of the primitive variables . To this end, let us consider a vector field that varies approximately linearly with distance as near . Up to linear order, the first derivative of is simply . The volumeaverage of the spatial derivatives of in the vicinity of is
(42) 
where we have assumed that the linear approximation is valid up to all the neighboring mesh generating points . It is straightforward to verify that the average matrix can be written as
(43) 
Writing the vector product in tensor form (where is a square matrix and and are vectors of dimension ), it is easy to prove the identity . Equivalently, going back to vector notation, we have , where, for simplicity, we used vector notation to denote a “cross product” between a vector and a matrix.
Therefore, the second term on the right hand side of Eq. (43) can be written as
(44) 
Here, the second term on the right hand side vanishes identically, because
(45) 
On the other hand, the first term on the right hand side of Eq. (44) can be rewritten by means of the replacement . Finally, identifying the vector with the gradient of a scalar quantity , and the matrix with the cellaveraged Hessian matrix , Eq. (44) takes the form
(46) 
The most noteworthy characteristic of this expression is that it is purely algebraic and explicit in nature. That is, the Hessian matrix of is simply a linear combination of the neighboring gradients in which the coefficients are predetermined quantities that depend only on the local mesh geometry. Each one of those neighboring gradients is, at the same time, a linear combination of its immediate neighbors’ scalar quantities (see Eq. 21 of Springel, 2010). Therefore, the Hessian estimate of Eq. (46) is a weighted linear combination of scalars from its immediate neighbors and from its neighbors’ neighbors and, as such, it implicitly employs a larger stencil than the one used for the gradients.
3.4 SlopeLimiting the Hessians
It is well known that higherorder reconstruction schemes are prone to produce spurious oscillations in the vicinity of steep gradients, unless this is prevented by appropriate slope limiter methods (Toro, 2009). These nonlinear corrections in the data reconstruction procedure prevent overshoots and allow for true discontinuities in the solutions. So far, we have discussed how to estimate second derivatives from first derivatives, which in turn are already reconstruction estimates obtained from the scalar physical quantities. Potential irregularities in the second derivative fields can lead to spurious oscillations and unphysical values of the viscous stress tensor at the cell boundaries. To alleviate this problem, we enforce local monotonicity of each component of the gradients, which is equivalent to smoothing out the Hessian estimates. In practice, this is achieved by replacing the Hessian matrix by a ‘slope limited’ version
(47) 
with
(48) 
are the slope limiters for each direction , and . This MINMODtype slopelimiting method is readily applicable for irregular meshes. The quantities are defined as
(49) 
where the are the components of the vector , i.e. the estimated change in the gradient between the centroid of the cell and the center of cell . The quantities and are the maximum and minimum of the th component of the cellcentered gradient estimates among all neighboring cells of cell , including itself.
To our knowledge, the slopelimiting technique has not been applied to the second derivatives before. However, its purpose is equivalent to the “flattening” procedure near shocks carried out by Colella & Woodward (1984) for parabolic reconstruction. In our approach, the suppression of oscillations near shocks is exclusively handled by the limitation of the gradients, since the reconstruction of hydrodynamic quantities is only of linear order. Thus, the Hessian limitation procedure serves the sole purpose of guaranteeing a smooth variation of the gradients and avoiding spuriously large viscous fluxes. Future improvements of the present method could however also employ these second derivatives for higherorder reconstructions of the scalars.
3.5 Time integration and timestep criterion
Because of the more complex mathematical properties of the NS equations compared with the Euler equations, obtaining a rigorous analytic expression analogous to the CFL stability criterion for the allowed time step size is not possible. However, MacCormack & Baldwin (1975) obtained an approximate semiempirical stability criterion when advective, viscous and heat diffusion terms are considered. When there is no heat flux, the timestep criterion can be written as (e.g. Kundu & Cohen, 2008)
(50) 
where is the standard CFLcriterion timestep except for the CourantFriedrichsLevy coefficient, which is absorbed into a “safety factor” (usually ). In Eq. (50), the cell Reynolds number is
(51) 
where is the velocity of the gas relative to the motion of the grid and is the effective radius of the cell, calculated as from the volume of a cell (or as from the area in 2D). Similar approaches to derive an appropriate NS timestep have also been described by Mavriplis & Jameson (1990) and Coirier & Powell (1996).
The numerical integration scheme we employ is time unsplit, that is, advective and diffusive fluxes are applied simultaneously during each hydrodynamic timestep and not sequentially (Eq. 25). The prediction stage, on the other hand, is operatorsplit, since the advective and diffusive terms are computed almost independently of each other. This is in part due to the nature of the standard onedimensional Riemann problem, whose solutions – strictly speaking – are only valid for the hyperbolic Euler problem, but are not solutions to the full NS equations with their additional parabolic terms. The validity of this approach ultimately relies on the assumption that the viscous terms in the NS equations are typically small perturbations to the Euler equations.
4 Numerical Test Results
To test the performance of AREPO when our new treatment of viscous fluxes is included, we have carried out a number of test simulations for physical situations with known analytic or quantitative solutions. Usually, the problems with known exact solutions are either of selfsimilar type or have symmetries that make the nonlinear term proportional to vanish identically. Owing to these limitations, numerical simulations of situations with experimentally wellestablished behavior, such as flow past a circular cylinder, have become commonplace in testing the performance of NS codes. We will therefore also carry out such qualitative benchmarks, besides looking at a few simple problems with analytic solutions.
4.1 Diffusion of a Vortex Sheet
A simple problem of laminar flow in the presence of viscosity is given by the vortex sheet diffusion test. In this problem, the initial velocity field at is given by with for and for . Because of the symmetry of the problem, the NS equations reduce to a 1D diffusion equation
(52) 
with solution (e.g. Kundu & Cohen, 2008)
(53) 
In Figure 5, we show the time evolution we obtain for a twodimensional simulation domain with initially uniform pressure and density (), and with a velocity field given by . The mesh generating points were distributed regularly at the initial time to produce a Cartesian mesh. As the system evolves, the velocity and the vorticity fields as a function of time and vertical coordinate follow the exact solution remarkably well. It is worth pointing out that the initial singularity in the vorticity field is unresolved numerically (and thus appears as being uniformly zero throughout the domain), since the system is started with an exact sharp discontinuity. Static, perfectly aligned meshes with slope limitation techniques will typically maintain this unresolved vorticity and thus no diffusion will proceed unless some numerical perturbations are seeded that break the mesh alignment of the initial state (a common way to overcome this difficulty is to start the system according to Eq. (53) at such that there is initial vorticity). However, the moving mesh of AREPO “sees” a nonzero velocity gradient as soon as the upper and lower halves of the domain become unaligned with respect to each other. This happens because, as soon as a cell shifts its position, the number of its neighbors that have a drastically different velocity increases and so does the “statistical weight” of the discontinuity. At this point, the slopelimiting technique, which had ignored the discontinuity in the perfectly aligned mesh, now identifies the local variation as “real” and the vorticity field is “detected”.
Fig. 6 shows the corresponding twodimensional velocity field of the diffusing vortex sheet test at four different times, together with the geometry of the underlying Voronoi mesh. The mesh geometry nicely shows how the cells transform from a Cartesian configuration to an unstructured mesh, while the velocity field evolves from a piecewise constant state with a central discontinuity to a smoothly varying shear flow due to the effects of viscosity.
4.2 Diffusion of a Gaussian Vortex
The twodimensional circular velocity distribution corresponding to an irrotational vortex of circulation is
(54) 
where the vorticity is zero everywhere except at the origin (, i.e. a vortex line). In a viscous fluid, this velocity profile has to be sustained by a point source of vorticity at the origin (e.g. an infinitely thin rotating cylinder) otherwise the vortex line will decay in a similar way as the vortex sheet in the previous example. If the velocity at the origin is set impulsively to zero, the subsequent evolution of the azimuthal velocity is given by
(55) 
while the vorticity evolves as
(56) 
and the Laplacian of the velocity field is
(57) 
Because of its geometry, this problem is significantly more challenging than the vortex sheet test considered above and cannot be impulsively started at precisely . Besides the initial singularity in the vorticity field, the velocity field is divergent as we approach the origin. In addition, it is not possible to capture the azimuthal velocity field when the distance from the origin is comparable to the grid resolution. At the same time, the azimuthal velocity field is challenging for the boundary conditions, because the problem is selfsimilar in nature and therefore natural boundaries do not exist. These problems did not exist for the vortexsheet problem, which is of onedimensional nature. Nevertheless, evolving the system from an initial time minimizes most of these complications. In addition, we extend the computational domain far beyond the region of interest, such that boundaries become essentially irrelevant during the timespan of the numerical solution.
We setup a Cartesian mesh () with an imposed initial velocity profile of
(58) 
corresponding to a Gaussian vortex that we center in the middle of the domain, which extends over the range , and thus accommodates a radial range from to . The adopted physical parameters are , , , and the initial density field is constant with . The pressure field, however, is not uniform because the fluid is not started from rest. We obtain the correct pressure profile from the radial component of the equation of motion:
and thus the initial pressure profile is
where is an integration constant. The precise value of is irrelevant for the similarity solution presented here, because it is obtained for incompressible flow. In our numerical experiments (which are compressible), we set such that at .
Fig. 7 shows the time evolution of the velocity field, the vorticity field and the Laplacian field for a Gaussian vortex started on an initially Cartesian mesh. We find not only that the velocity evolves as expected based on the similarity solution, but the first and second derivatives also show excellent agreement with the analytic expectations. These results validate both the space and timeaccuracy of our viscous integration scheme, as well as the accuracy with which the second derivatives are estimated.
4.3 Plane Poiseuille and Couette Flows
Next, we consider impulsivelystarted plane Poiseuille and Couette flows where a fluid between two parallel plates is initially at rest, and then, suddenly, either pressure gradients or plate motions are applied. The timedependent solution has the form , where the horizontal velocity can be decomposed into steady and timedependent parts, . In the presence of a pressure gradient and an upper plate moving at constant speed , the steady state solution is the wellknown expression
(59) 
for which the special cases and are commonly known as plane Poiseuille flow and plane Couette flow, respectively.