Viscous flow on a moving mesh

Multi-dimensional, compressible viscous flow on a moving Voronoi mesh

D. J. Muñoz V. Springel, R. Marcus, M. Vogelsberger and L. Hernquist
Harvard Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138
Heidelberg Institute for Theoretical Studies, Schloss-Wolfsbrunnenweg 35, 69118 Heidelberg, Germany
Zentrum für Astronomie der Universität Heidelberg, ARI, Mönchhofstr. 12-14, 69120 Heidelberg, Germany

Numerous formulations of finite volume schemes for the Euler and Navier-Stokes 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 mesh-generating 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 mesh-twisting effects are avoided. Our implementation considers the full Navier-Stokes 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 Navier-Stokes 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.

hydrodynamics – methods: numerical.
pagerange: LABEL:firstpageLABEL:lastpagepubyear: 2012

1 Introduction

The last two decades have seen remarkable advances in the numerical solution of the compressible Navier-Stokes (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 mesh-free 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, finite-volume implementations of the two-dimensional 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 mesh-generating 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 steady-state 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 well-defined 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, two-dimensional 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 three-dimensional 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 quasi-Lagrangian 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 mesh-generating 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 Voronoi-ALE 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 Voronoi-based finite volume scheme to the two-dimensional 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 self-gravity 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 Galilean-invariance 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 high-speed flows are present. While conventional fixed-mesh 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 Galilean-invariant solutions (modulo floating point round-off 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 magneto-hydrodynamic turbulence or anisotropic plasma viscosity. Global simulations of cold accretion disks around protostellar objects (e.g. see de Val-Borro et al., 2006) still include shear viscosity coefficients in the form of a Shakura-Sunyaev 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 Spitzer-Braginskii 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 well-specified 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 second-order 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 Navier-Stokes Equations

The compact form of the Euler equations, when written in terms of the vector of conserved quantities (Toro, 2009) is




and where


is the mass-momentum-energy 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


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,


needs to explicitly account for the work done by viscous forces.

A general parametrization of the viscous stress tensor is given by


Often, the viscous stress tensor is decomposed into a traceless part and a diagonal part, such that the first corresponds to constant-volume shear deformations (often called the rate-of-deformation 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 so-called 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 hard-sphere, Chapman-Enskog approach to kinetic theory, a non-zero value for the bulk viscosity is obtained. In an extension of the hard sphere fluid model, the Longuet-Higgins-Pople 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


where is a viscous source term given by


The solution of the Euler equations with source terms is often handled by operator-splitting methods (e.g. Toro, 2009; LeVeque, 2002). That is, the numerical scheme alternates between an advective step that solves the homogeneous part, and a source-term step. Thus, the solution of Eq. (7) is split into a two stage problem:


Typically, the source terms are more easily written in the primitive variable formulation of the Euler equations. A common choice of the primitive-variable 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 source-term step. The three-dimensional Euler equations can be written in the primitive variable form as (Toro, 2009)


For this choice of variables, the coefficient matrices are given by (Toro, 2009)


which is exactly equivalent to the familiar equations


In this formulation, the viscous terms of the NS equations, which affect only the velocity, are (e.g. Landau & Lifshitz, 1959)


An alternative to expressing the viscosity effects as source terms is to absorb them directly into the flux divergence,


which highlights the still conservative character of the NS equations. Here diffusive fluxes, defined by


are responsible for the effects of viscosity. An implementation of the diffusive fluxes in this conservation-law 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 MUSCL-Hancock 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 finite-volume form are

with (19)

where, in general, the intercell fluxes contain both advective and diffusive contributions,


The scheme used by AREPO is the finite volume MUSCL-Hancock approach, consisting of a MUSCL (Monotone Upstream-centered Schemes for Conservation Laws) linear reconstruction stage, and a Hancock two-stage time integration


where the numerical fluxes represent appropriately time-averaged approximations to the true flux across the interface shared by cells and . The time label in Eq. (21) indicates that an intermediate-stage (a half time-step evolution) has been performed to obtain the numerical estimate of , meaning that the time-stepping in Eq. (21) uses time-centered fluxes, giving it second-order accuracy. The Hancock part of the scheme is a two-step approach (the familiar predictor-corrector algorithm) in which the correction half-step is obtained from the solution of the 1-D Riemann problem across each face of the control volume. The general finite volume MUSCL-Hancock 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


where we denote by the estimated vector of conserved variables at the centroid of the -interface, obtained by linearly extrapolating the cell-centered values of the -th cell (on the “left” side) from , the cell’s center position, to . Similarly, corresponds to the estimates of the face-centroid values obtained by linear extrapolation of the cell-centered 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 time-step.

(II) Evolution of Boundary Extrapolated Values

This is, strictly speaking, the “predictor” half time-step. The conserved variables are evolved for with flux estimates obtained from the values at the beginning of the time-step:

(III) Solution of 1-D Riemann Problems and Computation of Godunov Fluxes

This corresponds to the “corrector” half time-step in the two-stage Hancock approach. Once the values to the right and left of the interface at time are known, the discontinuity is treated as a one-dimensional 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 time-centered fluxes used to update the system from the beginning of the time-step to its end,


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 Galilean-invariant discretization scheme is obtained in which the truncation error does not depend on the bulk velocity of the system.

Figure 1: Schematic representation of the mesh geometry and the MUSCL-Hancock integration scheme implemented in AREPO: a) TheVoronoi mesh is uniquely determined by the location of the mesh-generating points. b) A gradient estimate for all primitive variables is obtained from the immediate neighbors of a given cell. c) The gradient-estimation process is repeated for each cell in the domain and thus a piece-wise linear reconstruction is obtained for each primitive variable. d) The primitive variables are extrapolated toward each interface and evolved for half a time-step. e) For each face, a pair of extrapolated quantities for two neighboring cells and forms a local Riemann problem. f) The Riemann problem is solved for each face of a cell, yielding time-centered Godunov fluxes for the entire boundary of the control volume of cell . These fluxes are used for updating the conserved quantities of the cell through Eq. (21).
Figure 2: Detailed description of the flux calculation with a Riemann solver in step e) of Fig. 1. For the case of a moving mesh, the standard MUSCL-Hancock method needs to be augmented with Galilean-boosts, as described by Springel (2010): (1) The extrapolation towards each interface is followed by a Galilean boost of the velocities to the rest frame of the face, and by a rotation of the coordinate axes. Each face is then treated as a one-dimensional discontinuity. Thus, the axes are oriented in the rotated frame such that the -axis coincides with the normal to the face (left panel). (2) The primitive variables in the moving frame are evolved for half a time-step, including source terms if present (e.g. gravity or viscosity). (3) A one dimensional Riemann problem is solved at the interface. (4) The velocities are translated back to the lab frame and the advective fluxes are computed.

3.2 A MUSCL-Hancock Finite-Volume Scheme with Viscous Terms

A cell-centered, finite-volume solution of the NS equation can be written as


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 time-centered 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 time-centered fluxes, obtained here with the two-step Hancock technique, as described above. Thus, for estimating both and a half time-step predictor stage is required. In the MUSCL-Hancock approach for inviscid flow, this step is carried out by linear reconstruction from each cell center to the interface, followed by solving a one-dimensional 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 MUSCL-Hancock 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 cell-centered gradients linearly and evolve them for half a time-step.

  • Average the extrapolated velocity gradients at the interface and use them to estimate viscous fluxes.

To extrapolate the gradients from their cell-centered values to the interfaces, information about the higher-order 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 cell-centered 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 piece-wise 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 cell-centered values really represents the value at the midpoint of the two mesh-generating points, which, for a Voronoi mesh, can be substantially offset from the mid-point 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 three-stage 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 time-step evolution of the conserved quantities is carried out as in Eq. (25). However, the approach (A-C) is preferable to the Loh (2007) scheme because it uses time-centered 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 three-stage 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 MUSCL-Hancock 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 time-step evolution stage:


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 cell-centered source terms in Eq. (26) calls for second order derivatives of the velocity field. Thus, in addition to the cell-centered velocity gradients , and , the cell-centered 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 MUSCL-Hancock 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 slope-limited, 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 cell-centered estimates at each side of the face


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.

Figure 3: Schematic representation of the double linear reconstruction proposed in this work compared to standard linear reconstruction and parabolic reconstruction.

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 :


where is the Hessian matrix of the scalar quantity and is a tensor containing the third-order derivatives of (i.e. ). Truncating both Taylor expansions to first-order 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 cell-averaged over the mesh, the reconstruction procedure recovers the same polynomial. This condition is trivially satisfied for a linear reconstruction of the form . However, higher-order reconstruction schemes require the use of zero-mean polynomials, which, beyond first-order, 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 K-exact second-order 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 second-order-accurate) 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


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 time-step to obtain a time integration scheme that is consistent with the second-order accurate two-stage MUSCL-Hancock approach. In the latter, to extrapolate and evolve a scalar quantity we consider


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:


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):


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


where we introduced the rank- tensor . Since is a function of the primitive variables , the tensor can also be written as (see Appendix)


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 time-step evolved gradients of the velocity are (in analogy to Eq. 31):


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 time-centered 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.

Figure 4: Sketch illustrating the individual steps involved in the extrapolation and half time-step evolution of the gradients, analogous to the advective flux calculation shown in Fig. 2. The different steps are: (1) spatial extrapolation of the gradients, followed by (2) a time advance by according to Eq. (31), and (3) an approximate evaluation right at the interface. In step (4), the viscous fluxes are determined by evaluating Eq. (18) with the values of the primitive variables and the velocity gradients at the interface.

(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


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


where is the conserved variable vector at the centroid of the interface, as seen in the lab frame, obtained from the solution of a 1-D 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


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:


Just like with the advective fluxes, the flux tensor (Eq. 39) must be projected onto the normal of each -interface to obtain the 5-component vector


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 extrapolate-and-average 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 cell-centered 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 volume-average of the spatial derivatives of in the vicinity of is


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


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


Here, the second term on the right hand side vanishes identically, because


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 cell-averaged Hessian matrix , Eq. (44) takes the form


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 Slope-Limiting the Hessians

It is well known that higher-order 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 non-linear 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




are the slope limiters for each direction , and . This MINMOD-type slope-limiting method is readily applicable for irregular meshes. The quantities are defined as


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 cell-centered gradient estimates among all neighboring cells of cell , including itself.

To our knowledge, the slope-limiting 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 higher-order reconstructions of the scalars.

3.5 Time integration and time-step 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 semi-empirical stability criterion when advective, viscous and heat diffusion terms are considered. When there is no heat flux, the time-step criterion can be written as (e.g. Kundu & Cohen, 2008)


where is the standard CFL-criterion time-step except for the Courant-Friedrichs-Levy coefficient, which is absorbed into a “safety factor” (usually ). In Eq. (50), the cell Reynolds number is


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 time-step 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 time-step and not sequentially (Eq. 25). The prediction stage, on the other hand, is operator-split, since the advective and diffusive terms are computed almost independently of each other. This is in part due to the nature of the standard one-dimensional 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 self-similar type or have symmetries that make the non-linear term proportional to vanish identically. Owing to these limitations, numerical simulations of situations with experimentally well-established behavior, such as flow past a circular cylinder, have become common-place 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.

Figure 5: Diffusion of a vortex sheet. The two panels show the velocity along the -axis (left panel), and the vorticity (right panel), at times , 0.1, 0.2, 0.4, 0.8, and (from black to red), for a dynamic viscosity coefficient . The solid lines are given by the analytic solution described by Eqs. (53), while the solid circles are all 2500 cell-centered velocity and vorticity values of the initially Cartesian mesh. Note that the simulation is started with a sharp discontinuity in velocity and thus the -function vorticity field is initially unresolved. If the mesh would remain exactly Cartesian, the diffusion of vorticity would actually be suppressed in this case. Nevertheless, the small asymmetries introduced by the moving mesh trigger the diffusion regardless of the initially unresolved setup, and the time-dependent numerical result closely follows the expected exact solution.
Figure 6: Time evolution of the mesh geometry and the velocity field for a diffusing vortex sheet test. As the vorticity spreads from the center of the domain to the upper and lower boundaries, the mesh adapts to the continuous change in velocity until its original Cartesian structure disappears entirely. The color table (from blue to red) corresponds to the range between and in linear scale.

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


with solution (e.g. Kundu & Cohen, 2008)


In Figure 5, we show the time evolution we obtain for a two-dimensional 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 non-zero 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 slope-limiting 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 two-dimensional 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 piece-wise 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 two-dimensional circular velocity distribution corresponding to an irrotational vortex of circulation is


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


while the vorticity evolves as


and the Laplacian of the velocity field is


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 self-similar in nature and therefore natural boundaries do not exist. These problems did not exist for the vortex-sheet problem, which is of one-dimensional 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


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 time-accuracy of our viscous integration scheme, as well as the accuracy with which the second derivatives are estimated.

Figure 7: Time evolution of a diffusing Gaussian vortex. For each time (as labeled), we show the azimuthal velocity profile , the vorticity profile and the Laplacian profile , as computed by AREPO (blue points; only a random of the total shown) and compare it to the corresponding analytic expressions (solid red lines).
(a) Plane Poiseuille flow
(b) Plane Couette flow
Figure 8: Impulsively started plane Poiseuille and Couette flows as a function of time. a) Time evolution of the horizontal velocity profile versus vertical distance. Solid curves represent the analytic solutions of Eqs. (59) to (61) for and , at ten different times (time increasing from black to red). The data points correspond to all the cell-centered values of velocity along for a simulation started from rest. b) Time evolution of the horizontal velocity profile versus vertical distance. Solid curves represent the analytic solutions in Eqs. (59) to (61) for and at eight different times (time increasing from black to red). The data points correspond to all the cell-centered values of velocity along for a numerical simulation started from rest.
Figure 9: Time evolution of the mesh geometry and the velocity for flow between parallel plates. The horizontal velocity field for plane Poiseuille flow is rendered at four different times. The evolution of the velocity field (see Fig 8(a)) is accompanied by the evolution of the mesh from an initially Cartesian set up (top-left panel) to a fully unstructured grid by the time the flow has reached steady state (bottom-right panel). The (linear ) color scale ranges from blue () to red ().

4.3 Plane Poiseuille and Couette Flows

Next, we consider impulsively-started 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 time-dependent solution has the form , where the horizontal velocity can be decomposed into steady and time-dependent parts, . In the presence of a pressure gradient and an upper plate moving at constant speed , the steady state solution is the well-known expression


for which the special cases and are commonly known as plane Poiseuille flow and plane Couette flow, respectively.

The time dependent component is a solution of Eq. (52), subject to the initial condition and the boundary conditions at and . By separation of variables, the general solution is (e.g Graebel, 2007)


where the coefficients are determined by the initial condition