<figure><img src="image_1.png" /></figure>

HAL
open science

# Discontinuous Galerkin methods for axisymmetric flows

Anthony Bosco, Vincent Perrier

► To cite this version:

Anthony Bosco, Vincent Perrier. Discontinuous Galerkin methods for axisymmetric flows. Computers and Fluids, 2024, 270, pp.106139. 10.1016/j.compfluid.2023.106139. hal-04356856

HAL Id: hal-04356856

https://inria.hal.science/hal-04356856v1

Submitted on 20 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Discontinuous Galerkin methods for axisymmetric flows

Anthony Bosco, Vincent Perrier

October 13, 2023

## Abstract

In this article, high order discontinuous Galerkin methods for axisymmetric flows are developed. A first work is performed to put the equations in a so called canonical form to gather terms of similar nature, and also for ensuring the three dimensional conservativity. Then the numerical scheme is developed. Instead of relying on the strong axisymmetric formulation, we rely on the $\theta$-averaged form of a the standard three-dimensional discontinuous Galerkin method for the Navier-Stokes system. This allows to derive a numerical scheme without ambiguity on the source terms. The numerical scheme is new and makes link between different formulations previously proposed. A wide and progressive set of test cases, especially designed for axisymmetric problems is proposed. The numerical scheme is tested using these cases and the optimal order of accuracy is obtained.

# 1 Introduction

Discontinuous Galerkin methods allow for arbitrarily high order numerical approximations, but at the cost of increasing the number of degrees of freedom, in particular for three dimensional problems where the computational cost can become large for complicated problems, even when running on many CPUs and other high performance computing architectures.

When the geometries involved possess cylindrical symmetries, it is possible to cast the original three dimensional problem into a two dimensional problem, which then results in less expensive numerical simulations, which helps when considering the aforementioned high costs linked with high order approximations of DG methods.

Axisymmetric simplifications have been developed for the Euler equations in the context of finite volume method in [12, 24] and discontinuous Galerkin methods in [16], and the context of the discretization of compressible Navier-Stokes equations and RANS equations in [26] using the finite volume method, and in [3] in the context of discontinuous Galerkin methods. In these works, the numerical method is directly constructed from a strong form of the axisymmetric Navier-Stokes equations. Still, the derivation of a numerical method for such a strong form may be hazardous because the axisymmetric set of equations does not have a canonical form, and includes source terms that are not straightforward to integrate.

In [15, 22, 8, 25], a finite element discretization for the Laplace equations, Stokes equations and incompressible Navier Stokes equations is obtained by simplifying the weak formulation for a three dimensional numerical scheme using the axisymmetric assumption with cylindrical coordinates. A study of the numerical method obtained in these works concludes that the error estimates can be compared to classical finite element error estimates, but where the error is computed in a weighted norm, which appears naturally due to the metric $r$ in the integrals [9].

In this article, we are interested in extending the method of [15, 22, 8] to a discontinuous Galerkin method for compressible Navier-Stokes equations. The discontinuous Galerkin method
used is the BR2 discontinuous Galerkin method[4, 5]. The numerical scheme found is different from [3] and can be seen as an extension of [15, 22, 8, 25] to non-linear problems using the BR2 discontinuous Galerkin method. Also, as a source term depending on a gradient appears in the axisymmetric formulation, a link with the so-called dual consistent integration of source terms, which was introduced in [27] for RANS turbulence modeling, is performed. The last aim of this article is to provide a progressive set of test cases, including new ones, for testing the accuracy of the numerical method.

This article is organized as follows. In section 2 we start by a work on the different possible formulations of the compressible Navier-Stokes equations. We define canonical terms that can be clearly identified either by considering the conservativity of the equations or by highlighting the terms that cannot be formulated with so-called Cartesian operators.

Then in section 3, the numerical scheme for axisymmetric problems is derived by simplifying the three dimensional scheme written on a rotated two dimensional mesh. The numerical scheme obtained is compared with some previously proposed discretization of the canonical terms defined in section 2. The main difference from previous work is that this method of deriving the numerical scheme leads to the introduction of a weighted lifting operator. It also allows us to clarify the discretization of the canonical source term, as well as the discretization of the additional non-Cartesian flux, in particular in regards to the numerical fluxes in the BR2 discretization.

Last, in section 4, the numerical method obtained is validated using a range of test cases: scalar diffusion, scalar advection-diffusion [7], a laminar Poiseuille flow [14, 11] and a laminar flow past a sphere in the axisymmetric regime [13, 29, 17, 2].

# 2 On the different formulations of the equations

## 2.1 Cylindrical coordinates

Throughout this work, we consider geometries which possess rotational symmetries around an axis, which can be described using cylindrical coordinates. Let us consider the Cartesian basis ($e_x, e_y, e_z$) in a three dimensional space. Any vector $u$ can be decomposed into that basis

$$u = u_x e_x + u_y e_y + u_z e_z.$$

In order to represent the axisymmetric problem, we consider cylindrical coordinates in the moving basis ($e_r, e_\theta, e_z$), in which the coordinates ($r, \theta, z$) are defined, with reference to the Cartesian basis, with the relation

$$\begin{cases} x = r \cos \theta, \\ y = r \sin \theta, \\ z = z. \end{cases} \qquad (1)$$

Any vector in the Cartesian basis can be rewritten in cylindrical coordinates as

$$u = u_r e_r + u_\theta e_\theta + u_z e_z,$$

where the relation between the component of the vector are given by

$$\begin{aligned} u_x &= \cos \theta u_r - \sin \theta u_\theta \\ u_y &= \sin \theta u_r + \cos \theta u_\theta \\ u_z &= u_z, \end{aligned} \qquad (2)$$
or, introducing the rotation matrix R:

$$R(\theta) := \begin{pmatrix} \cos \theta & -\sin \theta & 0 \\ \sin \theta & \cos \theta & 0 \\ 0 & 0 & 1 \end{pmatrix}, \qquad (3)$$

the transformation between Cartesian and cylindrical coordinates is rewritten

$$\mathbf{u}^{\text{cart}} = R(\theta)\mathbf{u}^{\text{cyl}}.$$

Similarly, for any tensor $\bar{\mathbf{f}}$, the transformation is

$$\bar{\mathbf{f}}^{\text{cart}} = R(\theta)\bar{\mathbf{f}}^{\text{cyl}}R(\theta)^T.$$

Based on these transformations, the following definition of axisymmetric problems can be given **Definition 1** (Axisymmetry). A scalar is said to be axisymmetric if its value depends only on r and z.

A vector or a tensor are said to be axisymmetric if all their components in the cylindrical basis depend only on r and z.

A system of PDE is axisymmetric if all the scalars, vectors and tensors of the PDE are axisymmetric.

## 2.2 Cartesian formulation of axisymmetric Navier-Stokes equations

The Navier-Stokes equations reads, in three dimension

$$\left\{ \begin{array}{l} \partial_t \rho + \text{div}_x (\rho \mathbf{v}) = 0 \\ \partial_t (\rho \mathbf{v}) + \text{div}_x (\rho \mathbf{v} \mathbf{v}^T + P I_d) = \text{div}_x \bar{\tau} \\ \partial_t (\rho E) + \text{div}_x ((\rho E + P) \mathbf{v}) = \text{div}_x (-\mathbf{q} + \bar{\tau} \cdot \mathbf{v}), \end{array} \right. $$

where $\rho$ is the density, $\mathbf{v}$ the velocity, $P$ the pressure, $\bar{\tau}$ is the stress strain tensor:

$$\bar{\tau} = \mu (\nabla \mathbf{v} + \nabla^T \mathbf{v}) + \lambda (\text{div}_x \mathbf{v}) I_d,$$

$E$ is the specific total energy

$$E = \frac{\mathbf{v}^2}{2} + \varepsilon,$$

where $\varepsilon$ is the specific internal energy. Last, $\mathbf{q} = -\kappa \nabla T$ where $T$ is the temperature. The thermodynamic variables $\rho$, $P$, $T$ and $\varepsilon$ are linked with an equation of state: each thermodynamic variable can be expressed as a combination of two other independent thermodynamic variables. The Navier-Stokes equations will be shortly rewritten as

$$\partial_t \mathbf{U} + \text{div}_x \mathbf{F}_{\text{adv}}(\mathbf{U}) = \text{div}_x \mathbf{F}_{\text{diff}}(\mathbf{U}, \nabla \mathbf{U}),$$

where $\mathbf{F}_{\text{adv}}$ and $\mathbf{F}_{\text{diff}}$ are the advective and diffusive part of the Navier-Stokes equations

$$\mathbf{F}_{\text{adv}} = \begin{pmatrix} \rho \mathbf{v} \\ \rho \mathbf{v} \mathbf{v}^T + P I_d \\ (\rho E + P) \mathbf{v} \end{pmatrix} \qquad \mathbf{F}_{\text{diff}} = \begin{pmatrix} 0 \\ \bar{\tau} \\ \bar{\tau} \cdot \mathbf{v} - \mathbf{q} \end{pmatrix}.$$

As the diffusive flux is linear with respect to the gradient of $\mathbf{U}$, $\mathbf{F}_{\text{diff}}$ will be rewritten $\mathbf{A}(\mathbf{U})\nabla \mathbf{U}$ when needed.
Usually, an axisymmetric code is derived from a two dimensional Cartesian code. It may therefore be interesting to write the axisymmetric equations only in terms of Cartesian operators in the variables (r, z) in order to explicit the additional terms that appear in the axisymmetric case. For this, we first define the Cartesian operators as follows.

**Definition 2** (Cartesian operators). In the case with swirl, namely when $u_{\theta} \neq 0$, the following Cartesian operators in coordinates (r, z) are defined as

$$ 
\begin{aligned}
\nabla_{\{r,z\}}^{C} u &:= \left( \frac{\partial u}{\partial r}, 0, \frac{\partial u}{\partial z} \right)^{T} \\
\operatorname{div}_{\{r,z\}}^{C} \mathbf{u} &:= \frac{\partial \mathbf{u}_r}{\partial r} + \frac{\partial \mathbf{u}_z}{\partial z} \\
\operatorname{div}_{\{r,z\}}^{C} \bar{\mathbf{u}} &:= \left( \frac{\partial \bar{\mathbf{u}}_{rr}}{\partial r} + \frac{\partial \bar{\mathbf{u}}_{rz}}{\partial z} \right) e_r + \left( \frac{\partial \bar{\mathbf{u}}_{\theta r}}{\partial r} + \frac{\partial \bar{\mathbf{u}}_{\theta z}}{\partial z} \right) e_{\theta} + \left( \frac{\partial \bar{\mathbf{u}}_{rz}}{\partial r} + \frac{\partial \bar{\mathbf{u}}_{zz}}{\partial z} \right) e_z \\
\nabla_{\{r,z\}}^{C} \mathbf{u} &:= \begin{pmatrix} \frac{\partial \mathbf{u}_r}{\partial r} & \frac{\partial \mathbf{u}_{\theta}}{\partial r} & \frac{\partial \mathbf{u}_z}{\partial r} \\ 0 & 0 & 0 \\ \frac{\partial \mathbf{u}_r}{\partial z} & \frac{\partial \mathbf{u}_{\theta}}{\partial z} & \frac{\partial \mathbf{u}_z}{\partial z} \end{pmatrix}.
\end{aligned}
$$

In the case without swirl, namely when $u_{\theta} = 0$, the following Cartesian operators in coordinates (r, z) are defined as

$$ 
\begin{aligned}
\nabla_{\{r,z\}}^{C} u &:= \left( \frac{\partial u}{\partial r}, \frac{\partial u}{\partial z} \right)^{T} \\
\operatorname{div}_{\{r,z\}}^{C} \mathbf{u} &:= \frac{\partial \mathbf{u}_r}{\partial r} + \frac{\partial \mathbf{u}_z}{\partial z} \\
\operatorname{div}_{\{r,z\}}^{C} \bar{\mathbf{u}} &:= \left( \frac{\partial \bar{\mathbf{u}}_{rr}}{\partial r} + \frac{\partial \bar{\mathbf{u}}_{rz}}{\partial z} \right) e_r + \left( \frac{\partial \bar{\mathbf{u}}_{rz}}{\partial r} + \frac{\partial \bar{\mathbf{u}}_{zz}}{\partial z} \right) e_z \\
\nabla_{\{r,z\}}^{C} \mathbf{u} &:= \begin{pmatrix} \frac{\partial \mathbf{u}_r}{\partial r} & \frac{\partial \mathbf{u}_z}{\partial r} \\ \frac{\partial \mathbf{u}_r}{\partial z} & \frac{\partial \mathbf{u}_z}{\partial z} \end{pmatrix}.
\end{aligned}
$$

Considerations on the three dimensional conservativity of equations (see [10]) lead to the following natural definition

**Definition 3** (Canonical axisymmetric conservative form). A scalar conservation law

$$ \partial_t u + \operatorname{div}_x \mathbf{f} = 0 $$

has the following Canonical axisymmetric conservative form

$$ \partial_t (ru) + \operatorname{div}_{\{r,z\}}^{C} (rf) = 0. $$

A vectorial conservation law

$$ \partial_t u + \operatorname{div}_x \bar{\mathbf{f}} = 0, $$

has the following Canonical axisymmetric conservative form

$$ \partial_t (ru) + \operatorname{div}_{\{r,z\}}^{C} (rf) = \bar{f}_{\theta\theta} e_r - \bar{f}_{r\theta} e_{\theta}. $$

Based on the **Definition 3**, the Navier-Stokes system can be put in the following canonical axisymmetric conservative form, which is proven in [10].
**Proposition 1** (Canonical Cartesian axisymmetric conservative form). The axisymmetric Navier-Stokes equations can be put in the following canonical Cartesian axisymmetric conservative form

$$
\left\{
\begin{array}{l}
\partial_t (r\rho) + \operatorname{div}_{\{r,z\}}^C (r\rho\mathbf{v}) = 0 \\
\partial_t (r\rho\mathbf{v}) + \operatorname{div}_{\{r,z\}}^C (r(\rho\mathbf{v}\mathbf{v}^T + P\mathrm{I}_d)) = \operatorname{div}_{\{r,z\}}^C (r\bar{\tau}^C) + S_{[\rho\mathbf{v}]}^{NC} + \mathbf{S}^{\text{Can}} \\
\partial_t (r\rho E) + \operatorname{div}_{\{r,z\}}^C (r(\rho E + P)\mathbf{v}) = \operatorname{div}_{\{r,z\}}^C (r\kappa\nabla_{\{r,z\}}^C T) + \operatorname{div}_{\{r,z\}}^C (r\bar{\tau}^C \cdot \mathbf{v}) + S_{[\rho E]}^{NC},
\end{array}
\right.
$$

• In the case with swirl, the source terms are given by

$$
\mathbf{S}^{\text{Can}} = \begin{pmatrix} P + \rho \mathbf{v}_{\theta}^2 - \frac{2\mu + \lambda}{r} \mathbf{v}_r - \lambda \operatorname{div}_{\{r,z\}}^{C} \mathbf{v} \\ \mu \left( \frac{\partial \mathbf{v}_{\theta}}{\partial r} - \frac{\mathbf{v}_{\theta}}{r} \right) - \rho \mathbf{v}_r \mathbf{v}_{\theta} \\ 0 \end{pmatrix}
$$
$$
S_{[\rho \mathbf{v}]}^{NC} = \operatorname{div}_{\{r,z\}}^{C} \bar{\tau}^{NC} = \begin{pmatrix} \frac{\partial (\lambda \mathbf{v}_r)}{\partial r} \\ -\frac{\partial (\mu \mathbf{v}_{\theta})}{\partial r} \\ \frac{\partial (\lambda \mathbf{v}_r)}{\partial z} \end{pmatrix}
$$
$$
S_{[\rho E]}^{NC} = \operatorname{div}_{\{r,z\}}^{C} (\bar{\tau}^{NC} \mathbf{v}) = \frac{\partial}{\partial r} (\lambda \mathbf{v}_r^2 - \mu \mathbf{v}_{\theta}^2) + \frac{\partial}{\partial z} (\lambda \mathbf{v}_r \mathbf{v}_z).
$$

• In the no swirl case, this reduces to

$$
\left\{
\begin{aligned}
\mathbf{S}^{\text{Can}} &= \begin{pmatrix} P - \frac{2\mu + \lambda}{r} \mathbf{v}_r - \lambda \operatorname{div}_{\{r,z\}}^{C} \mathbf{v} \\ 0 \end{pmatrix} \\
S_{[\rho \mathbf{v}]}^{NC} &= \operatorname{div}_{\{r,z\}}^{C} \bar{\tau}^{NC} = \operatorname{div}_{\{r,z\}}^{C} (\lambda \mathbf{v}_r \mathrm{I}_d) = \nabla_{\{r,z\}}^{C} (\lambda \mathbf{v}_r) \\
S_{[\rho E]}^{NC} &= \operatorname{div}_{\{r,z\}}^{C} (\bar{\tau}^{NC} \mathbf{v}) = \operatorname{div}_{\{r,z\}}^{C} (\lambda \mathbf{v}_r \mathbf{v})
\end{aligned}
\right.
$$

The equations in Proposition 1 can be shortly rewritten as

$$
\begin{equation}
\begin{split}
\partial_t (r\mathbf{U}) + \operatorname{div}_{\{r,z\}}^C (r\mathbf{F}_{\text{adv}}^C(\mathbf{U})) ={}& \operatorname{div}_{\{r,z\}}^C (r\mathbf{F}_{\text{diff}}^C(\mathbf{U}, \nabla_{\{r,z\}}^C \mathbf{U})) \\
& + \operatorname{div}_{\{r,z\}}^C (\mathbf{F}_{\text{diff}}^{NC}(\mathbf{U})) + \mathbf{S}^{\text{Can}}(\mathbf{U}, \nabla_{\{r,z\}}^C \mathbf{U}).
\end{split}
\tag{4}
\end{equation}
$$

We remark that:

• The first source term, $\mathbf{S}^{\text{Can}}$ is compatible with (3). This source term depends on the variables and their first order derivative.

* By using only two dimensional Cartesian operators in the $(r, z)$ variables, another source term induced by the second order terms appears, $\operatorname{div}_{\{r,z\}}^{C} \mathbf{F}_{\text{diff}}^{NC}(\mathbf{U})$. This source term depends also on the variables and their first order derivative.

These two source terms essentially come from different transformations of the equations. Even if they depend on the same variables, we will see that their discretization should be performed differently. Next section is dedicated to the discretization of these operators.

3 Derivation of the numerical scheme

3.1 Three dimensional numerical scheme

The aim of this section is to recall the discontinuous Galerkin numerical scheme for the Navier-
Stokes system. A possible discretization of the Navier-Stokes system is the so-called BR2
formulation which was first introduced in [6], and later in [4, 5], and analyzed in [1] for elliptic systems. The formulation of this numerical scheme is briefly recalled. For the sake of simplicity, we will not consider the boundary conditions for the moment.

The computational domain $\Omega$ is divided into a partition of cells $\hat{K} \in \mathcal{K}$. We denote by $\hat{\mathcal{F}}$ the set of faces, namely the $(d-1)$ dimensional intersections of elements of $\mathcal{K}$. A face linking the elements $\hat{K}_i$ and $\hat{K}_j$ is denoted by $\hat{\Gamma}_{ij}$. It is oriented such that $\hat{K}_i$ is on the left, and $\hat{K}_j$ is on the right of $\hat{\Gamma}_{ij}$. A set of basis $\hat{\Psi}$ is defined on each cell $\hat{K}$; the $\hat{\Psi}$ are continuous inside $\hat{K}$, and a priori discontinuous at the faces.

Each face $\mathcal{F}$ is oriented. From this orientation, the left cell (denoted with a $L$ subscript) is defined as the cell such that the normal to the face is outgoing, and the other one is defined as the right cell (denoted with a $R$ subscript). We then define the *average* operator on each face as

$$ \{\hat{\Psi}\} = \frac{\hat{\Psi}_L + \hat{\Psi}_R}{2}, $$

and the *jump* operator on each face as

$$ \llbracket \hat{\Psi} \rrbracket = \hat{\Psi}_L - \hat{\Psi}_R. $$

Then the BR2 formulation reads

$$ \sum_{\hat{K} \in \mathcal{K}} \left( \int_{\hat{K}} \hat{\Psi} \partial_t \mathbf{U} - \int_{\hat{K}} \mathbf{F}_{adv}(\mathbf{U}) \nabla \hat{\Psi} + \int_{\hat{K}} (\mathbf{A}(\mathbf{U}) \nabla \mathbf{U}) \nabla \hat{\Psi} \right) \\ + \sum_{\hat{\Gamma}_{ij} \in \mathcal{F}} \left( \int_{\hat{\Gamma}_{ij}} \tilde{F}(\mathbf{U}, \hat{n}) \llbracket \hat{\Psi} \rrbracket - \int_{\hat{\Gamma}_{ij}} \left( \{\mathbf{A}(\mathbf{U}) \nabla \mathbf{U}\} \llbracket \hat{\Psi} \rrbracket + \{\mathbf{A}^T(\mathbf{U}) \nabla \hat{\Psi}\} \llbracket \mathbf{U} \rrbracket \right) \right) \\ - \sum_{\hat{\Gamma}_{ij} \in \mathcal{F}} \left( \int_{\hat{\Gamma}_{ij}} \eta \{\mathbf{A}(\mathbf{U}) \ell^{ij} (\llbracket \mathbf{U} \rrbracket)\} \llbracket \hat{\Psi} \rrbracket \right) = 0. \quad (5) $$

This formulation depends on the definition of a *local lifting operator*, which is defined for a scalar $u$ as

$$ \forall \hat{\Psi} \quad \int_{\hat{K}_i \cup \hat{K}_j} \hat{\Psi} \cdot \hat{\ell}^{ij} (\llbracket u \rrbracket) = - \int_{\hat{\Gamma}_{ij}} \{\{\hat{\Psi}\}\} \cdot \llbracket u \hat{n} \rrbracket, \quad (6) $$

and which is defined componentwise by using (6) for a vector. Defining the *global* lifting operator $\mathcal{L}$ in each cell as

$$ \forall \hat{K}_i \in \mathcal{K} \quad \hat{\mathcal{L}}(\llbracket \mathbf{U} \rrbracket) = \sum_{j \in \mathcal{N}(i)} \hat{\ell}^{ij}(\llbracket \mathbf{U} \rrbracket), \quad (7) $$

the numerical scheme (5) can be rewritten as

$$ \sum_{\hat{K} \in \mathcal{K}} \left( \int_{\hat{K}} \hat{\Psi} \partial_t \mathbf{U} - \int_{\hat{K}} \mathbf{F}_{adv}(\mathbf{U}) \nabla \hat{\Psi} + \int_{\hat{K}} (\mathbf{A}(\mathbf{U}) (\nabla \mathbf{U} + \hat{\mathcal{L}}(\llbracket \mathbf{U} \rrbracket))) \nabla \hat{\Psi} \right) \\ + \sum_{\hat{\Gamma}_{ij} \in \mathcal{F}} \left( \int_{\hat{\Gamma}_{ij}} \tilde{F}(\mathbf{U}, \hat{n}) \llbracket \hat{\Psi} \rrbracket - \int_{\hat{\Gamma}_{ij}} \left( \{\mathbf{A}(\mathbf{U}) (\nabla \mathbf{U} + \eta \hat{\ell}^{ij}(\llbracket \mathbf{U} \rrbracket))\} \llbracket \hat{\Psi} \rrbracket \right) \right) = 0. \quad (8) $$

The last term in (5) and (8) depends on a penalty parameter $\eta$, which is set to ensure stability of discretization. Indeed, it can be shown (e.g.[1]) that the BR2 method is stable for

$$ \eta \geq N_{\partial}, \quad (9) $$

where $N_{\partial}$ is the maximum number of face of the cells adjacent to the face $\hat{\Gamma}_{ij}$.
<figure><img src="image_3.png" /><figcaption>Figure 1: Three dimensional mesh cell generated by the rotation of a triangle with an angle $\theta_0$.</figcaption></figure>

## 3.2 Derivation of the numerical scheme for axisymmetric models

As was remarked in the section 2, the axisymmetric model when written with Cartesian operators depends on source term that include derivatives, for which the derivation of an accurate numerical scheme is not straightforward [27]. Also the radial coordinate $r$ appears in front of the source term and the fluxes, which suggests that the system of equations is written on a space with a metric $r$ rather than on a Cartesian metric. This raises some questions on the expression of the lifting operator, that were already addressed for linear problems [15, 22, 9, 8, 25].

For avoiding any ambiguity in the derivation of the numerical scheme for axisymmetric models, the derivation of the numerical scheme will not rely directly on the continuous axisymmetric formulation of the Navier-Stokes system, but rather on the axisymmetric simplification of the three dimensional numerical scheme (5) or (8), with the lifting operator defined by (6), (7). More precisely, starting from a computational domain $\Omega$ in the $(r, z)$ plane, a three dimensional computational domain $\hat{\Omega}$ is generated by rotation around the $z$ axis. We also start with a mesh of $\Omega$, composed of cells $K \in \mathcal{K}$. A mesh of $\hat{\Omega}$ is obtained by rotating the cells $K \in \mathcal{K}$, see Figure 1. We denote by $\mathcal{F}$ the set of the faces of $\mathcal{K}$, by $\Gamma_{ij}$ the face linking $K_i$ and $K_j$. A discontinuous approximation space is defined on the mesh $\mathcal{K}$: continuous in each cell $K \in \mathcal{K}$, and discontinuous across the faces $\Gamma_{ij}$. A set of basis $\hat{\Psi}$ is defined on each cell $\hat{K}$, such that

$$ \forall (r, \theta, z) \in \hat{K} \quad \hat{\Psi}(r, \theta, z) = \Psi(r, z). $$

Then (5) or (8) provide a suitable numerical scheme for the three dimensional equations. Suppose now that no variable is depending on $\theta$. Then all the integrals on the cells $\hat{K}_i \in \hat{K}$ can be simplified into integrals on $K_i \in \mathcal{K}$, and all the integrals on $\hat{\Gamma}_{ij}$ can be simplified into integrals on $\Gamma_{ij}$. The remaining part of this section is dedicated to the reduction of all the integrals of (5) or (8) into integrals on the faces and cells of $\mathcal{K}$.

### 3.2.1 Simplification of the mass integral

The mass integral is

$$ \int_{\hat{K}} \hat{\Psi} \partial_t U d\hat{K}. $$

It can be integrated in cylindrical coordinates. Provided $U$ depends only on $r$ and $z$ (namely, the vectors were expressed in the basis $e_r, e_\theta$), we find

$$ \int_{\hat{K}} \hat{\Psi} \partial_t U d\hat{K} = 2\pi \int_{K} \Psi \partial_t U r dK. $$
### 3.2.2 Simplification of the lifting operator

On each face, the local lifting operator is defined as

$$ \forall \hat{\Psi} \quad \int_{\hat{K}_i \cup \hat{K}_j} \hat{\Psi} \cdot \hat{\ell}^{ij} (\llbracket u \rrbracket) = - \int_{\hat{\Gamma}_{ij}} \{\{\hat{\Psi}\}\} \cdot \llbracket u \hat{n} \rrbracket. $$

$\hat{\Psi}$ does not depend on $\theta$. As $\hat{n}$ is orthogonal to $e_\theta$, the right hand side always vanishes when $\hat{\Psi}$ is parallel to $e_\theta$. This means that $\hat{\ell}^{ij}(\llbracket u \rrbracket)$ has no component on $e_\theta$. Supposing that $u$ is axisymmetric, everything can be integrated in cylindrical coordinates for finding

$$ \forall \hat{\Psi} \quad 2\pi \int_{K_i \cup K_j} \hat{\Psi} \cdot \hat{\ell}^{ij} (\llbracket u \rrbracket) r \, dK = -2\pi \int_{\Gamma_{ij}} \{\{\hat{\Psi}\}\} \cdot \llbracket u \hat{n} \rrbracket r \, d\Gamma. $$

This last equality leads to the following definition

**Definition 4 (Weighted lifting operator).** The local weighted lifting operator $\mathcal{L}_w^{ij}(\llbracket u \rrbracket)$ is a three dimensional vector defined as

$$ \mathcal{L}_w^{ij}(\llbracket u \rrbracket) := (\mathcal{L}_{w,r}^{ij}(\llbracket u \rrbracket), 0, \mathcal{L}_{w,z}^{ij}(\llbracket u \rrbracket))^T, $$

with $\mathcal{L}_{w,k}^{ij}(\llbracket u \rrbracket)$ defined for $k=r$ or $k=z$ as

$$ \forall \Psi \quad \int_{K_i \cup K_j} \Psi \mathcal{L}_{w,k}^{ij}(\llbracket u \rrbracket) r \, dK = - \int_{\Gamma_{ij}} \{\{\Psi\}\} \llbracket u \rrbracket \mathbf{n} \cdot \mathbf{e}_k r \, d\Gamma. $$

The global weighted lifting operator $\mathcal{L}_w(\llbracket u \rrbracket)$ is a three dimensional vector defined as

$$ \mathcal{L}_w(\llbracket u \rrbracket) = \sum_{j \in \mathcal{N}(i)} \mathcal{L}_w^{ij}(\llbracket u \rrbracket). $$

Note that for a vector $\mathbf{u}$, the weighted lifting tensor is a 3 × 3 matrix defined as

$$ \mathcal{L}_w^{ij}(\llbracket \mathbf{u} \rrbracket) := (\mathcal{L}_w^{ij}(\llbracket \mathbf{u}_r \rrbracket), \mathcal{L}_w^{ij}(\llbracket \mathbf{u}_\theta \rrbracket), \mathcal{L}_w^{ij}(\llbracket \mathbf{u}_z \rrbracket)). $$

**Remark 1.** The weighted lifting operator defined here is derived from the original three dimensional BR2 formulation. We note that, in previous work on axisymmetric problems for the compressible Navier-Stokes equations [3], the lifting operator used in the BR2 method is the usual Cartesian operator. Note that the numerical scheme with the weighted lifting operator is directly inheriting all the usual properties of the three dimensional scheme (consistency, stability, coercivity).

### 3.2.3 Simplification of the cell integrals

**Proposition 2 (Axisymmetric form of the cell integrals).** The cell integrals in the three dimensional BR2 formulation

$$ \frac{1}{2\pi} \int_{\hat{K}} (\mathbf{F}(\mathbf{U}) - \mathbf{A}(\mathbf{U}) (\nabla \mathbf{U} + \hat{\mathcal{L}}(\llbracket \mathbf{U} \rrbracket))) \nabla \hat{\Psi} \, d\hat{K} $$

can be simplified into

$$ \int_K r (\mathbf{F}_{adv}^C(\mathbf{U}) - \mathbf{F}_{diff}^C(\mathbf{U}, \nabla_{\{r,z\}}^C \mathbf{U} + \mathcal{L}_w(\llbracket \mathbf{U} \rrbracket))) \nabla_{\{r,z\}}^C \Psi \, dK \\ - \int_K \mathbf{F}_{diff}^{NC}(\mathbf{U}) \nabla_{\{r,z\}}^C \Psi \, dK - \int_K \mathbf{S}^{Can}(\mathbf{U}, \nabla_{\{r,z\}}^C \mathbf{U} + \mathcal{L}_w(\llbracket \mathbf{U} \rrbracket))) \Psi \, dK. \quad (10) $$
Proof. **Continuous operator** We first deal with the integral

$$ \int_K (\mathbf{F}(\mathbf{U}) - \mathbf{A}(\mathbf{U})\nabla\mathbf{U}) \nabla\hat{\Psi}. $$

The mass component expressed in cylindrical coordinates gives

$$ (\rho \mathbf{v}) \cdot \nabla \hat{\Psi} = \rho \mathbf{v}_r \frac{\partial \Psi}{\partial r} + \rho \mathbf{v}_z \frac{\partial \Psi}{\partial z} = \rho \mathbf{v} \cdot \nabla_{\{r,z\}}^{C} \Psi. $$

Concerning the momentum component, we have

$$ \begin{aligned} (\rho v \mathbf{v}^T + P I_d - \bar{\tau}) : \nabla \hat{\Psi} &= (\rho v \mathbf{v}^T + P I_d - \bar{\tau}) : \begin{pmatrix} \frac{\partial \hat{\Psi}_r}{\partial r} & \frac{\partial \hat{\Psi}_\theta}{\partial r} & \frac{\partial \hat{\Psi}_z}{\partial r} \\ -\frac{\partial \hat{\Psi}_\theta}{\partial r} & \frac{\partial \hat{\Psi}_r}{\partial r} & 0 \\ \frac{\partial \hat{\Psi}_r}{\partial z} & \frac{\partial \hat{\Psi}_\theta}{\partial z} & \frac{\partial \hat{\Psi}_z}{\partial z} \end{pmatrix} \\ &= (\rho v \mathbf{v}^T + P I_d - \bar{\tau}) : \nabla_{\{r,z\}}^C \hat{\Psi} - \frac{\hat{\Psi}_\theta}{r} (\rho v_r \mathbf{v}_\theta - \bar{\tau}_{r\theta}) \\ &\quad + \frac{\hat{\Psi}_r}{r} (P + \rho v_\theta^2 - \bar{\tau}_{\theta\theta}) \\ &= (\rho v \mathbf{v}^T + P I_d - \bar{\tau}) : \nabla_{\{r,z\}}^C \hat{\Psi} + \frac{1}{r} \mathbf{S}^{\text{Can}} \cdot \hat{\Psi} \\ (\rho v \mathbf{v}^T + P I_d - \bar{\tau}) : \nabla \hat{\Psi} &= (\rho v \mathbf{v}^T + P I_d - \bar{\tau}^C) : \nabla_{\{r,z\}}^C \hat{\Psi} - \frac{1}{r} \bar{\tau}^{\text{NC}} : \nabla_{\{r,z\}}^C \hat{\Psi} + \frac{1}{r} \mathbf{S}^{\text{Can}} \cdot \hat{\Psi}. \end{aligned} $$

It remains to deal with the energy equation, for which we have

$$ \begin{aligned} &((\rho E + P) \mathbf{v} - \kappa \nabla T - \bar{\tau} \cdot \mathbf{v}) \cdot \nabla \hat{\Psi} \\ &= ((\rho E + P) - \kappa \nabla_{\{r,z\}}^C T - \bar{\tau}^C \cdot \mathbf{v}) \cdot \nabla_{\{r,z\}}^C \hat{\Psi} - \frac{1}{r} (\bar{\tau}^{\text{NC}} \cdot \mathbf{v}) \cdot \nabla_{\{r,z\}}^C \hat{\Psi}. \end{aligned} $$

We can now integrate everything in cylindrical coordinates, replace the three dimensional test functions (with hats) by two dimensional test functions for finding

$$ \frac{1}{2\pi} \int_{\hat{K}} (\mathbf{F}(\mathbf{U}) - \mathbf{A}(\mathbf{U})\nabla\mathbf{U}) \nabla\hat{\Psi} = \int_{K} r (\mathbf{F}_{adv}^{C}(\mathbf{U}) - \mathbf{F}_{diff}^{C}(\mathbf{U}, \nabla_{\{r,z\}}^{C}\mathbf{U})) \nabla_{\{r,z\}}^{C} \Psi \, dK \\ - \int_{K} \mathbf{F}_{diff}^{NC}(\mathbf{U}) \nabla_{\{r,z\}}^{C} \Psi \, dK - \int_{K} \mathbf{S}^{\text{Can}}(\mathbf{U}, \nabla_{\{r,z\}}^{C}\mathbf{U}) \Psi \, dK. $$

**Lifting part** We are now interested in simplifying the following cell integral

$$ \int_{\hat{K}} (\mathbf{A}(\mathbf{U})\hat{\mathcal{L}}(\llbracket \mathbf{U} \rrbracket)) \nabla \hat{\Psi}. $$

The mass equation has no component. Concerning the momentum equation, we need to deal with

$$ \bar{\tau}(\llbracket \mathbf{u} \rrbracket) := \mu(\hat{\mathcal{L}}(\llbracket \mathbf{u} \rrbracket) + \hat{\mathcal{L}}^T(\llbracket \mathbf{u} \rrbracket)) + \lambda \mathrm{Tr}(\hat{\mathcal{L}}(\llbracket \mathbf{u} \rrbracket)). $$

It is worth noting that as the second row of $\bar{\tau}(\llbracket \mathbf{u} \rrbracket)$ is zero, we have $\bar{\tau}(\llbracket \mathbf{u} \rrbracket) = \bar{\tau}^C(\llbracket \mathbf{u} \rrbracket)$. Then

$$ \begin{aligned} \bar{\tau}(\llbracket \mathbf{u} \rrbracket) : \nabla \hat{\Psi} &= \bar{\tau}(\llbracket \mathbf{u} \rrbracket) : \begin{pmatrix} \frac{\partial \hat{\Psi}_r}{\partial r} & \frac{\partial \hat{\Psi}_\theta}{\partial r} & \frac{\partial \hat{\Psi}_z}{\partial r} \\ -\frac{\partial \hat{\Psi}_\theta}{\partial r} & \frac{\partial \hat{\Psi}_r}{\partial r} & 0 \\ \frac{\partial \hat{\Psi}_r}{\partial z} & \frac{\partial \hat{\Psi}_\theta}{\partial z} & \frac{\partial \hat{\Psi}_z}{\partial z} \end{pmatrix} \\ &= \bar{\tau}^C(\llbracket \mathbf{u} \rrbracket) : \nabla_{\{r,z\}}^C \hat{\Psi} - \frac{\hat{\Psi}_\theta}{r} \bar{\tau}_{r\theta}(\llbracket \mathbf{u} \rrbracket) + \frac{\Psi_r}{r} \bar{\tau}_{\theta\theta}(\llbracket \mathbf{u} \rrbracket) \end{aligned} \quad (11) $$
In this last equality, we recognize the lifting of the Cartesian diffusive operator, and the lifting evaluated at the components of the Canonical source term that depend on derivatives.

The energy term can be dealt in the same manner, and the Cartesian diffusive operator evaluated with the lifting is found. As done in the continuous case, the integrals can then be evaluated in cylindrical coordinates, and be explicitly integrated in $\theta$.

Last, the lifting operator is replaced by the weighted lifting operator of Definition 4.

**Summary** Gathering the part that depends on continuous derivatives and the part that depends on the lifting operator, we finally find that the three dimensional integral can be reduced into the computation of several two dimensional integrals which is given in (10) □

### 3.2.4 Simplification of the face integrals

For simplifying the face integrals, we will require the following property on the hyperbolic numerical flux

**Definition 5** (Rotation invariant numerical flux). A numerical flux $\tilde{F}(U_L, U_R, \hat{n})$ for the Euler system is invariant by rotation if for any $R$, rotation of angle $\theta$, and by denoting

$$ \mathcal{R} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & R & 0 \\ 0 & 0 & 1 \end{pmatrix}, $$

the following holds for any rotation

$$ \mathcal{R}^T \tilde{F}(\mathcal{R}U_L, \mathcal{R}U_R, R\hat{n}) = \tilde{F}(U_L, U_R, \hat{n}). $$

**Example** Let us consider the Lax-Friedrich numerical flux

$$ \tilde{F}(U, \hat{n}) = \{\mathbf{F}_{adv}(U) \cdot \hat{n}\} + \lambda \llbracket \mathbf{U} \rrbracket \quad (12) $$

At the interface, using the definition of the normal vector, the normal flux $\mathbf{F} \cdot \mathbf{n}$, is given, in cylindrical coordinates, by

$$ \mathbf{F}_{adv}(\mathbf{U}) \cdot \hat{\mathbf{n}} = \begin{pmatrix} \rho u_n \\ \rho u_n u_r + P n_r \\ \rho u_n u_\theta \\ \rho u_n u_z + P n_z \\ (\rho E + P) u_n \end{pmatrix}, $$

where $u_n = u_r n_r + u_z n_z$ given the definition of the normal vector $\hat{n}$. Equation (12) in axisym- metric form is then given by

$$ \tilde{F}(\mathbf{U}, \hat{\mathbf{n}}) = \left\{ \left\{ \begin{pmatrix} \rho u_n \\ \rho u_n u_r + P n_r \\ \rho u_n u_\theta \\ \rho u_n u_z + P n_z \\ (\rho E + P) u_n \end{pmatrix} \right\} \right\} + \lambda \llbracket \mathbf{U} \rrbracket. $$

Note that, in the case without swirl, the third component is removed and gives the classical 2D Lax-Friedrich numerical flux.

We are now able to show the following proposition.
**Proposition 3** (Axisymmetric form of the face integrals). If the numerical flux is invariant by rotation, then the face integrals in the three dimensional BR2 formulation

$$ \frac{1}{2\pi} \left( \int_{\hat{\Gamma}_{ij}} \tilde{F}(\mathbf{U}, \hat{\mathbf{n}}) [\hat{\mathbf{\Psi}}] - \int_{\hat{\Gamma}_{ij}} \left( \left\{ \mathbf{A}(\mathbf{U}) (\nabla \mathbf{U} + \eta \hat{\mathbf{\ell}}^{ij} ([\mathbf{U}])) \right\} \right) [\hat{\mathbf{\Psi}}] \right), $$

can be simplified into the axisymmetric form

$$ \int_{\Gamma_{ij}} \tilde{F}(\mathbf{U}, \hat{\mathbf{n}}) [\hat{\mathbf{\Psi}}] \, \mathrm{rd}\Gamma_{ij} - \int_{\Gamma_{ij}} \left\{ \mathbf{F}_{\mathrm{diff}}^{C}(\mathbf{U}, \nabla \mathbf{u} + \eta \hat{\mathbf{\ell}}_{w}^{ij} ([\mathbf{U}])) \right\} [\hat{\mathbf{\Psi}} \otimes \mathbf{n}] \, \mathrm{rd}\Gamma_{ij} \\ - \int_{\Gamma_{ij}} \left\{ \mathbf{F}^{NC}(\mathbf{U}) \right\} [\hat{\mathbf{\Psi}} \otimes \mathbf{n}] \, \mathrm{d}\Gamma_{ij}. \qquad (13) $$

*Proof.* The face $\hat{\Gamma}_{ij}$ in the three dimensional formulation is induced by the rotation of a face $\Gamma_{ij}$, such that, using cylindrical coordinates, we may write

$$ d\hat{\Gamma}_{ij} = r d\theta d\Gamma_{ij}, $$

where $d\Gamma_{ij}$ is the infinitesimal surface element of the face in the two dimensional mesh. Moreover, on a rotated mesh cell, the normal vector $\hat{\mathbf{n}}$ to a face $\hat{\Gamma}_{ij}$ does not have an angular component $n_{\theta}$ in cylindrical coordinates, i.e.

$$ \hat{\mathbf{n}} = (n_r, 0, n_z)^T. $$

**Convective part** First, let us consider the convective part, i.e.

$$ \int_{\hat{\Gamma}_{ij}} \tilde{F}(\mathbf{U}, \hat{\mathbf{n}}) [\hat{\mathbf{\Psi}}]. $$

Considering a rotated mesh cell, the face integral for the convective part can be rewritten as

$$ \int_{\hat{\Gamma}_{ij}} \tilde{F}(\mathbf{U}, \hat{\mathbf{n}}) [\hat{\mathbf{\Psi}}] = \int_{0}^{2\pi} \int_{\Gamma_{ij}} \tilde{F}(\mathbf{U}, \hat{\mathbf{n}}) [\Psi] \, \mathrm{rd}\theta \, \mathrm{d}\Gamma_{ij}. \qquad (14) $$

With some abuse of notation, we denote by $\tilde{F}$ both the numerical flux that use the whole solution $\mathbf{U}: \tilde{F}(\mathbf{U}, \hat{\mathbf{n}})$ and the one that depends on the left and right states $\tilde{F}(\mathbf{U}_L, \mathbf{U}_R, \hat{\mathbf{n}})$.

In the special case of Definition 5 where the rotation is done around the z-axis, this means that if a rotation invariant numerical flux is used, then its value in cylindrical coordinates does not depend on $\theta$ provided $\mathbf{U}_L$ and $\mathbf{U}_R$ do not depend on $\theta$ in cylindrical coordinates. Then (14) can be simplified as

$$ \int_{\hat{\Gamma}_{ij}} \tilde{F}(\mathbf{U}, \hat{\mathbf{n}}) [\hat{\mathbf{\Psi}}] = \int_{0}^{2\pi} \int_{\Gamma_{ij}} \tilde{F}(\mathbf{U}, \hat{\mathbf{n}}) [\Psi] \, \mathrm{rd}\theta \, \mathrm{d}\Gamma_{ij} \qquad (15) \\ = 2\pi \int_{\Gamma_{ij}} \tilde{F}(\mathbf{U}, \hat{\mathbf{n}}) [\Psi] \, \mathrm{rd}\Gamma_{ij}. $$

**Diffusive part** Let us now consider the diffusive part of the BR2 formulation, which, in three dimensions, can be written

$$ \int_{\hat{\Gamma}_{ij}} \left\{ \mathbf{F}^v (\mathbf{U}, \nabla \mathbf{U} + \eta \hat{\mathbf{\ell}}^{ij} ([\mathbf{U}])) \right\} [\hat{\mathbf{\Psi}} \otimes \mathbf{n}]. $$
Similarly as in the case of cell integrals, the viscous flux can be written in cylindrical coordinates as a combination of a two dimensional conservative part, and a "non conservative" part, i.e.

$$F^v(\mathbf{U}, \nabla \mathbf{U} + \eta \hat{\mathbf{l}}^{ij}(\llbracket \mathbf{U} \rrbracket)) = F_C^v(\mathbf{U}, \nabla \mathbf{U} + \eta \hat{\mathbf{l}}^{ij}(\llbracket \mathbf{U} \rrbracket)) + F_{NC}^v(\mathbf{U}),$$

where the second part only depends on $\mathbf{U}$, and not on the gradient and lifting operator.

The mass equation has no component for this term. We are interested in the momentum equation, in which the viscous flux can be rewritten as

$$\tau(\mathbf{U}, \nabla \mathbf{U} + \eta \hat{\mathbf{l}}^{ij}(\llbracket \mathbf{U} \rrbracket)) = \tau^C(\mathbf{U}, \nabla \mathbf{U} + \eta \hat{\mathbf{l}}^{ij}(\llbracket \mathbf{U} \rrbracket)) + \frac{1}{r} \tau^{NC}(\mathbf{U}),$$

Direct computation of $\hat{\Psi} \otimes \mathbf{n}$ in a rotated mesh cell, in cylindrical coordinates, gives

$$\hat{\Psi} \otimes \mathbf{n} = \begin{pmatrix} \Psi_r n_r & 0 & \Psi_r n_z \\ \Psi_\theta n_r & 0 & \Psi_\theta n_z \\ \Psi_z n_r & 0 & \Psi_z n_z \end{pmatrix} = \Psi \otimes \mathbf{n},$$

since the angular component $n_\theta = 0$.

Using cylindrical coordinates, we can then simplify the integral in a similar way as in the convective flux, i.e.

$$\int_{\hat{\Gamma}_{ij}} \{\hat{\tau}\} : [\hat{\Psi} \otimes \mathbf{n}] d\hat{\Gamma}_{ij} = \int_{0}^{2\pi} \int_{\Gamma_{ij}} \{\tau\} : [\Psi \otimes \mathbf{n}] r d\theta d\Gamma_{ij} \qquad (16) \\ = 2\pi \int_{\Gamma_{ij}} \{\tau\} : [\Psi \otimes \mathbf{n}] r d\Gamma_{ij}.$$

Similarly as in the case of cell integrals, we define

$$\eta \bar{\tau}(\llbracket \mathbf{u} \rrbracket) := \eta \left( \mu (\hat{\ell}(\llbracket \mathbf{u} \rrbracket) + \hat{\ell}^T(\llbracket \mathbf{u} \rrbracket)) + \lambda \mathrm{Tr}(\hat{\ell}(\llbracket \mathbf{u} \rrbracket)) \right),$$

such that the lifting part $\eta \tau(\llbracket \mathbf{u} \rrbracket)$ can be written in a similar way as previous face integrals, i.e.

$$\int_{\hat{\Gamma}_{ij}} \{\eta \bar{\tau}^c(\llbracket \mathbf{u} \rrbracket)\} : [\hat{\Psi} \otimes \mathbf{n}] d\hat{\Gamma}_{ij} = \int_{0}^{2\pi} \int_{\Gamma_{ij}} \{\eta \bar{\tau}^c(\llbracket \mathbf{u} \rrbracket)\} : [\Psi \otimes \mathbf{n}] r d\theta d\Gamma_{ij} \qquad (17) \\ = 2\pi \int_{\Gamma_{ij}} \{\eta \bar{\tau}^c(\llbracket \mathbf{u} \rrbracket)\} : [\Psi \otimes \mathbf{n}] r d\Gamma_{ij}.$$

The non-Cartesian part of the stress tensor $\tau^{NC}$ can be treated similarly, resulting, for the momentum part, in the following

$$\frac{1}{2\pi} \int_{\hat{\Gamma}_{ij}} \{\hat{\tau}\} : [\hat{\Psi} \otimes \mathbf{n}] = \int_{\Gamma_{ij}} \{\tau^C\} : [\Psi \otimes \mathbf{n}] r d\Gamma_{ij} + \int_{\Gamma_{ij}} \frac{1}{r} \{\tau^{NC}\} : [\Psi \otimes \mathbf{n}] r d\Gamma_{ij}. \qquad (18)$$

Regarding the energy equation, proceeding similarly as in the case of the cell integrals, we compute

$$\int_{\hat{\Gamma}_{ij}} \{(-\kappa \nabla T - \bar{\tau} \cdot \mathbf{v}) \cdot \mathbf{n}\} [\hat{\Psi}] d\hat{\Gamma}_{ij} = 2\pi \int_{\Gamma_{ij}} \{(-\kappa \nabla^C_{\{r,z\}} T - \tau^c \cdot \mathbf{v}) \cdot \mathbf{n}\} [\Psi] r d\Gamma_{ij} \qquad (19) \\ -2\pi \int_{\Gamma_{ij}} \frac{1}{r} \{\tau^{NC} \cdot \mathbf{v} \} \cdot \mathbf{n} [\Psi] r d\Gamma_{ij}.$$

**Final expression for the face integrals** Gathering (15), (16), (17), (18), and (19), the face integrals in the original three dimensional formulation can be rewritten as (13)

□
**Remark 2.** The product $\hat{\Psi} \otimes n$ can be rewritten

$$ \begin{pmatrix} \hat{\Psi}_r n_r & 0 & \hat{\Psi}_r n_z \\ \hat{\Psi}_\theta n_r & 0 & \hat{\Psi}_\theta n_z \\ \hat{\Psi}_z n_r & 0 & \hat{\Psi}_z n_z \end{pmatrix} = \underbrace{\begin{pmatrix} \hat{\Psi}_r n_r & 0 & \hat{\Psi}_r n_z \\ 0 & 0 & 0 \\ \hat{\Psi}_z n_r & 0 & \hat{\Psi}_z n_z \end{pmatrix}}_{2d \text{ part}} + \underbrace{\begin{pmatrix} 0 & 0 & 0 \\ \hat{\Psi}_\theta n_r & 0 & \hat{\Psi}_\theta n_z \\ 0 & 0 & 0 \end{pmatrix}}_{\text{swirl part}} $$

Then, using this decomposition, the product $\tau : (\hat{\Psi} \otimes n)$ can be rewritten as a two dimensional $(r, z)$ part and a swirl contribution

$$ \tau : (\hat{\Psi} \otimes n) = (\tau : \Psi \otimes n)_{\{r,z\}} + (\tau : \Psi \otimes n)_\theta, $$

in accordance with the strong form of the axisymmetric equations (4).

### 3.2.5 Final expression of the scheme

Gathering the results of the previous subsections, namely the equations (10), (13), the following proposition was proven

**Proposition 4** (Axisymmetric discontinuous Galerkin scheme). If a two dimensional axisymmetric model is discretized on a three dimensional mesh obtained by rotation of a two dimensional mesh, then the final numerical scheme on the two dimensional mesh reads

$$ \begin{aligned} & \int_K \Psi \partial_t U \,rdK + \int_K (\mathbf{F}_{adv}^C(\mathbf{U}) - \mathbf{F}_{diff}^C(\mathbf{U}, \nabla_{\{r,z\}}^C \mathbf{U} + \mathcal{L}_w(\llbracket \mathbf{U} \rrbracket))) \nabla_{\{r,z\}}^C \Psi \,rdK \\ & \quad - \int_K \mathbf{F}_{diff}^{NC}(\mathbf{U}) \nabla_{\{r,z\}}^C \Psi \,dK - \int_{\Gamma_{ij}} \tilde{F}(\mathbf{U}, \hat{\mathbf{n}}) \llbracket \Psi \rrbracket \,rd\Gamma_{ij} \\ & \quad + \int_{\Gamma_{ij}} \{\{\mathbf{F}_{diff}^C(\mathbf{U}, \nabla \mathbf{U} + \eta \boldsymbol{\ell}_w^{ij}(\llbracket \mathbf{U} \rrbracket))\}\} \llbracket \Psi \otimes \mathbf{n} \rrbracket \,rd\Gamma_{ij} + \int_{\Gamma_{ij}} \{\{\mathbf{F}^{NC}(\mathbf{U})\}\} \llbracket \Psi \otimes \mathbf{n} \rrbracket \,d\Gamma_{ij} \\ & = \int_K \mathbf{S}^{Can}(\mathbf{U}, \nabla_{\{r,z\}}^C \mathbf{U} + \mathcal{L}_w(\llbracket \mathbf{U} \rrbracket)) \Psi \,rdK. \end{aligned} \quad (20) $$

## 3.3 Remark on the discretization of the right hand side in Proposition 1

The DG discretization obtained by simplifying the original three dimensional problem into a two dimensional problem using rotated mesh cells results in a numerical scheme which is comparable to that obtained by directly obtaining a weak formulation from the strong form of the axisymmetric Navier Stokes equations (4).

However, in the axisymmetric Navier-Stokes equations, the integration of the non Cartesian part of the diffusive flux is not straightforward, since it is a first order term, and a numerical flux has to be specified in the face integrals,

$$ \int_{\Gamma_{ij}} \hat{\mathbf{F}}^{NC}(\mathbf{U}) \llbracket \Psi \otimes \mathbf{n} \rrbracket \,d\Gamma. $$

The numerical scheme obtained by simplifying the three dimensional problem directly introduces a centered numerical flux

$$ \int_{\Gamma_{ij}} \{\{\mathbf{F}^{NC}(\mathbf{U})\}\} \llbracket \Psi \otimes \mathbf{n} \rrbracket \,d\Gamma. $$
It can also be noticed that the source term $S^{Can}$ is composed of two parts. The first part is the direct integration on cells of the source term appearing in the axisymmetric equations,

$$ \int_K S^{Can}(U, \nabla_{\{r,z\}} U) \Psi r dK. $$

We notice however that an additional term also appears, which depends of the lifting operator and takes the form

$$ \int_K S^{Can}(0, \mathcal{L}_w(\llbracket \mathbf{U} \rrbracket)) \Psi r dK. $$

The relevance of the lifting operator in the source term will be discussed in the next section.

## 3.4 Link with dual consistent and asymptotic dual consistent discretization of the source term

Dual consistency is a key property in the analysis of discontinuous Galerkin methods [1], for instance in the case of elliptic problem, it is used to prove optimal order of convergence in $L_2$ norm, and dual consistency of various DG methods have been studied in [19].

In the context of RANS turbulence modeling, Oliver [27] extended this analysis to the source terms arising from turbulence models, which depends non linearly on the gradient of the solution $\nabla u$. In particular, using a simplified elliptic model, he showed that in order to obtain a discretization which converges with optimal order of convergence, a simple integration of the source term

$$ \int_K S(U, \nabla U) \Psi dK, $$

does not result in a dual consistent discretization, and that in order to recover that property, correction terms need to be added. In particular, he showed that adding the lifting operator from the BR2 discretization to the source term

$$ \int_K S(U, \nabla U + \mathcal{L}(\llbracket \mathbf{U} \rrbracket)) \Psi dK, $$

results in an asymptotically dual consistent method, and recovers optimal order of convergence in $L_2$ norm.

The discretization obtained when simplifying the original three dimensional formulation can be further developed. Indeed, the source term, which is composed of two parts - one from the convective flux, and one from the diffusive flux - can be rewritten in a more compact way. Indeed, the diffusive part of the source term is given by

$$ \int_K S^v(U, \nabla U) \Psi r dK = \int_K \frac{\hat{\Psi}_r}{r} (\bar{\tau}_{\theta\theta}(\nabla U) + \bar{\tau}_{\theta\theta}(\llbracket \mathbf{u} \rrbracket)) r dK, $$

where the first part depends on the gradient $\nabla U$ and the second part only on the lifting operator, as defined (11). Since the stress tensor is linear, we find that

$$ \int_K S^v(U, \nabla U) r dK = \int_K \frac{\hat{\Psi}_r}{r} (\bar{\tau}_{\theta\theta}(\nabla U + \mathcal{L}(\llbracket \mathbf{U} \rrbracket))) r dK, $$

which is the asymptotically dual consistent discretization from [27]. Thus, the simplification of the original three dimensional BR2 DG formulation using a rotated mesh is nothing but a dual consistent discretization of the canonical source term.
## 3.5 Boundary conditions

In discontinuous Galerkin methods, boundary conditions are imposed weakly through the boundary face integrals. For geometries in which the axis of symmetry is not part of the domain, for instance in the case of a flow between two cylinders, classical boundary conditions on the flow variables U can be used.

For geometries in which the axis of symmetry $r = 0$ is part of the computational domain, the boundary conditions at the axis need to be specified. In particular, it is readily observable that most boundary integrals vanish for faces in which $r = 0$. Only the non-Cartesian part of the stress tensor remain, i.e.

$$ \int_{\Gamma_{ij}} \{\{\mathbf{F}^{NC}(\mathbf{U})\}\} : [\Psi \otimes \mathbf{n}] d\Gamma. $$

$\mathbf{F}^{NC}(\mathbf{U})$ depends only on $v$. On the axis boundary, $v_{\theta}$ should be set to 0 whereas the values of $v_r$ and $v_z$ should be taken as the inner values.

## 3.6 Practical implementation

The implementation started from a code able to deal with one, two, and three dimensional Navier-Stokes equations discretized with the BR2 discontinuous Galerkin method on unstructured hybrid meshes in equation (5)

The code is organized in modules called cells and faces integrators able to compute the integrals of (5).

The derivation of the numerical scheme for the axisymmetric Navier-Stokes equations in (20) resulted in a formulation that is slightly different from the original Cartesian formulation. As a consequence, in order to implement axisymmetric models on top of the existing two dimensional Cartesian code, the following change were considered:

* Add the weight $r$ in each of the integrals to switch from two dimensional Cartesian models to axisymmetric models.

* Add the weight $r$ also in the integrals required for the computation of the lifting for changing the classical lifting operator into the weighted lifting operator.

* Adding axisymmetric correction terms in the momentum and energy equations in the cells and faces which were not originally present in the two dimensional Cartesian formulation but appeared when simplifying the three dimensional problem (namely new terms arising from the non-Cartesian flux), i.e.

$$ \int_K \mathbf{F}^{NC}(\mathbf{U}) \nabla \Psi \, dK, $$

for the cell integrals. They are added to the cell integrals in the BR2 formulation

$$ \sum_{K \in K} \left( \int_K \Psi \partial_t \mathbf{U} r dK - \int_K (\mathbf{F}_{adv}(\mathbf{U}) - (\mathbf{A}(\mathbf{U})\nabla \mathbf{U})) \nabla \Psi r dK + \int_K \mathbf{F}^{NC}(\mathbf{U}) \nabla \Psi dK \right) $$

Similarly, in the case of the face integrals, the derivation of the scheme starting from the three dimensional problem resulted in an additional non-Cartesian flux at the interface similar to the original diffusive flux, i.e. a centered flux

$$ \int_{\Gamma_{ij}} \{\{\mathbf{F}^{NC}(\mathbf{U})\}\} : [\Psi \otimes \mathbf{n}] d\Gamma_{ij}. $$
This term is added to the face integrals in the BR2 formulation

$$ \sum_{\Gamma_{ij} \in \mathcal{F}} \left( \int_{\Gamma_{ij}} (\tilde{F}(\mathbf{U}, \mathbf{n})[\Psi] - (\{\mathbf{A}(\mathbf{U})\nabla \mathbf{U}\}\Psi] + \{\mathbf{A}^T(\mathbf{U})\nabla \Psi\}[\mathbf{U}]) \right) r d\Gamma_{ij} + \int_{\Gamma_{ij}} \{\mathbf{F}^{\text{NC}}(\mathbf{U})\} : [\Psi] d\Gamma_{ij} $$

* For models in which the transformation from Cartesian to axisymmetric results in additional source terms, the source term cell integrals are computed in a dual consistent way by adding the global lifting operator to the gradient

$$ \int_K \mathbf{S}(\mathbf{U}, \nabla \mathbf{U} + \mathcal{L}_w(\llbracket \mathbf{U} \rrbracket)) \Psi r dK. $$

**Remark 3.** Using the formulation in (5), the global lifting operator needs not be explicitly computed and stored in memory. However, when needed, it can be reconstructed from the local lifting operators using the definition in equation (7), where the local lifting operators are computed in the penalty face integrals.

# 4 Numerical results

The validity and accuracy of the axisymmetric models implemented in the code AeroSol are assessed through a range of test cases of increasing difficulty: scalar diffusion, scalar advection-diffusion, Poiseuille flow, and a flow past a sphere at low Reynolds number. When an exact solution is known, the weighted $L_2$ error is computed, i.e.

$$ ||u - u_h||_r^2 = \int_{\Omega} |u - u_h|^2 r d\Omega. $$

This weighted error can be interpreted as the three dimensional error, up to a factor $2\pi$ on a rotated mesh. Moreover, in the case of a scalar elliptic problems, it can be seen that weighted function spaces appear as the natural space to study the numerical method [15, 22, 9, 8], and that typical error estimates can be obtained in weighted $H^1$ and $L_2$ norms. Considering the properties discussed in Remark 1, this weighted norm is the natural norm in which all the convergence theorems for linear problems hold.

In these simulations, the penalty coefficient used in the BR2 discretization is the one in (9). In particular, since quadrangular elements are used, we have $N_\partial = 4$ and the penalty coefficient is set to $\eta = 6$. In the case of piecewise constant approximations, the penalty term also acts as a consistency term. In that case, the penalty parameter is set to $\eta = 2$ for quadrangular cells.

## 4.1 Scalar diffusion

As a first validation test case, a scalar elliptic problem is considered, in which no additional source needs to be added when switching from a three dimensional problem to the axisymmetric model. For an axisymmetric geometry, the diffusion equations reads

$$ -\frac{1}{r} \mathrm{div}_{\{r,z\}}^{\mathcal{C}} (r \nabla_{\{r,z\}}^{\mathcal{C}} u) = f(r, z). $$
<figure><img src="image_4.png" /><figcaption>Figure 2: Convergence of the weighted $L_2$ in the case of an axisymmetric scalar diffusion problem</figcaption></figure>

We consider the axisymmetric diffusion problem in a unit square $\Omega: (r, z) \in [0, 1] \times [0, 1]$, and the right hand side is set to

$$f(r, z) = \frac{\sin(r)}{r} \exp(-z).$$

Dirichlet boundary conditions are chosen such that the exact solution takes the form

$$u(r, z) = \cos(r) \exp(-z).$$

The diffusion problem is solved using uniform quadrangular meshes with size $h_0$ with $h_0 = 0.1$ and $i = 0, 1, 2, 3$, i.e. with $10 \times 10, 20 \times 20, 40 \times 40$ and $80 \times 80$ cells in the $r$ and $z$ direction. The numerical results obtained for the weighted $L_2$ error are shown in Figure 2a in the case the weighted lifting operator is used, and in Figure 2b with the regular lifting operator. In both cases, the optimal order of accuracy $O(h^{k+1})$ of the numerical scheme for $k = 0, 1, 2, 3$ is obtained. In particular, from a numerical point of view, we observe that both lifting operators lead to very similar results. This is not surprising for $k \ge 1$ as the role of the lifting operator in (20) is to stabilize the numerical scheme and because the two formulations of the lifting operators are equivalent up to the multiplication or division by the radius of a cell or a face. Still, the weighted lifting operator is the only one up to our knowledge on which the choice of penalization parameter (9) gives a guarantee of optimal convergence for the axisymmetric formulation on linear problems. In the following, the weighted lifting operator, as determined when deriving the numerical scheme, will be used.

## 4.2 Advection-Diffusion problem

We consider the advection-diffusion problem in a hollow cylinder, adapted from the three dimensional test case in [7]. In axisymmetric form, the scalar advection-diffusion writes

$$\beta \cdot \nabla_{\{r,z\}}^{C}(ru) + \operatorname{div}_{\{r,z\}}^{C} \cdot (r\mu\nabla u) = 0,$$

where the advection velocity $\beta = (0, \beta_z)^T$. Dirichlet boundary conditions are imposed in order to obtain an exact solution of the form

$$u(r, z) = \frac{\log(r)}{\log(r_2)} \left[ k_1 \exp\left(\frac{\beta_z z}{\mu}\right) + k_2 \right],$$
<figure><img src="image_5.png" /><figcaption>Figure 3: Convergence of the weighted error with the mesh size h for the axisymmetric advection-diffusion problem.</figcaption></figure>

where $k_1$ and $k_2$ are two constants given by

$$ \left\{ \begin{aligned} k_1 &= \frac{1}{2 \left[ 1 - \exp \left( \frac{\beta_z L}{\mu} \right) \right]} \\ k_2 &= 1 - k_1, \end{aligned} \right. $$

where L is the height of the cylinder. Dirichlet boundary conditions are used in accordance with the exact solution $u(r, z)$.

The convergence of the axisymmetric advection-diffusion problem is analyzed in a rectangular domain of dimension $\Omega: (r, z) \in [0.5, 1] \times [0, 1]$ using uniformly refined quadrangular meshes with size $h_0/2^i$, with $i = 0, 1, 2, 3$, i.e. with $5 \times 10, 10 \times 20, 20 \times 40$ and $40 \times 80$ cells in the r and z direction. The parameters are set as follows: the axial advection velocity is set to $\beta_z = 20$ and the diffusion coefficient $\mu = 5$.

The weighted $L_2$ error as a function of the mesh size h for orders $k = 0, 1, 2, 3, 4$ is shown in Figure 3. The order of convergence is also computed, and the optimal order of convergence of the method $O(h^{k+1})$ is obtained.

## 4.3 Poiseuille flow

In order to validate the axisymmetric solver for the Navier Stokes equations without swirl, we consider the case of a Poiseuille flow in an axisymmetric domain, which can be separated into two configurations: a flow between two cylinders and a flow in a tube.

These test-cases were adapted from three dimensional cylindrical Poiseuille flows implemented in the library AeroSol, which were previously used as validation test-cases for high order meshes [14]. These test-cases are also based on the two dimensional Poiseuille flow problem introduced in [11].

In both cases, the flow is driven by an external body force $f$, which induces a right hand side

$$ \left( \begin{array}{c} f_{\rho} \\ f_{\rho u} \\ f_{\rho c} \end{array} \right) = \left( \begin{array}{c} 0 \\ \mathbf{f} \\ \mathbf{f} \cdot \mathbf{u} \end{array} \right), $$
where the term $f \cdot u$ is the work of the external force. The force is chosen so as to impose a parabolic velocity profile, i.e. the following equations are being solved

$$ \partial_t (r\mathbf{U}) + \mathrm{div}_{\{r,z\}}^C (r\mathbf{F}_{\mathrm{adv}}^C(\mathbf{U})) = \mathrm{div}_{\{r,z\}}^C (r\mathbf{F}_{\mathrm{diff}}^C(\mathbf{U}, \nabla_{\{r,z\}}^C \mathbf{U})) \\ + \mathrm{div}_{\{r,z\}}^C (\mathbf{F}_{\mathrm{diff}}^{NC}(\mathbf{U})) + \mathbf{S}^{\mathrm{Can}}(\mathbf{U}, \nabla_{\{r,z\}}^C \mathbf{U}) + r\rho \begin{pmatrix} f_{\rho} \\ \mathbf{f} \\ \mathbf{f} \cdot \mathbf{u} \end{pmatrix}. $$

### 4.3.1 Poiseuille flow between two cylinders

We first consider a flow between two cylinders of radii $r_1$ and $r_2$, in which the velocity $u(r, z) = u_z(r)e_z$ is assumed to be parabolic, i.e. of the form

$$ u_z(r) = u_0 \left(1 - \frac{r}{r_1}\right) \left(1 - \frac{r}{r_2}\right), $$

which can be obtained by balancing the axial momentum equation so that

$$ f_z = -\frac{\partial}{\partial r} \left( \mu r \frac{\partial u_z}{\partial r} (r) \right) $$

The forcing term $f$ can be computed and is given by

$$ \mathbf{f} = \begin{pmatrix} f_r \\ f_z \end{pmatrix} = \begin{pmatrix} 0 \\ -\frac{\mu u_0}{\rho} \frac{4r-r_1-r_2}{r_1 r_2 r} \end{pmatrix}, $$

such that the right hand side is given by

$$ \begin{pmatrix} f_{\rho} \\ \mathbf{f} \\ f_{\rho e} \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ -\frac{\mu u_0}{\rho} \frac{4r-r_1-r_2}{r_1 r_2 r} \mathbf{f} \cdot \mathbf{u} \end{pmatrix}. $$

The cylinders are considered as isothermal walls with fixed temperature $T_w$,

$$ T(r_1) = T(r_2) = T_w, $$

and periodic boundary conditions are used for in the streamwise direction.

Under these assumptions, the exact solution for the temperature is [14]

$$ T(r) = \frac{4 \mathrm{Pr} \mathrm{Ma}_{\infty}^2 (\gamma - 1) T_w}{9 (\Delta r)^4} \left[ 9 \left( (r_1 + r_2)^2 + r^2 \right) r^2 - 16 (r_1 + r_2) r^3 \right] + k_1 + k_2 \log(r). $$

with $\Delta r = r_2 - r_1$, and where $k_1$ and $k_2$ are two constants determined so as to impose the wall temperature $T_w$, with

$$ \left\{ \begin{array}{l} k_1 = \frac{(\log(r_2) - \log(r_1))T_w - \log(r_2)T_a + \log(r_1)T_b}{\log(r_2) - \log(r_1)}, \\ k_2 = \frac{T_a - T_b}{\log(r_2) - \log(r_1)}, \end{array} \right. $$

where the temperatures $T_a$ and $T_b$ are defined as

$$ \left\{ \begin{array}{l} T_a = \frac{4 \mathrm{Pr} \mathrm{Ma}_{\infty}^2 (\gamma - 1) T_w}{9 (\Delta r)^4} \left[ 9 \left( (r_1 + r_2)^2 + r_1^2 \right) r_1^2 - 16 (r_1 + r_2) r_1^3 \right], \\ T_b = \frac{4 \mathrm{Pr} \mathrm{Ma}_{\infty}^2 (\gamma - 1) T_w}{9 (\Delta r)^4} \left[ 9 \left( (r_1 + r_2)^2 + r_1^2 \right) r_2^2 - 16 (r_1 + r_2) r_2^3 \right]. \end{array} \right. $$
<figure><img src="image_6.png" /><figcaption>Figure 4: Exact solution for the Poiseuille flow between two cylinders.</figcaption></figure>

The parameters are set as follows: we consider a flow of air, treated as an ideal gas with heat capacity ratio $\gamma = 1.4$. The centerline velocity $u_c = 1$, the reference density $\rho_{ref} = 0.001$, the Reynolds number $Re = 1$, the wall temperature $T_w = 10$, the Mach number $Ma_\infty = 0.2$ and the Prandtl number $Pr = 0.7$. Since we consider the centerline velocity $u_c$ as a parameter, we may also rewrite the velocity

$$u_z(r) = -\frac{4u_c r_1 r_2}{(\Delta r)^2} \left(1 - \frac{r}{r_1}\right) \left(1 - \frac{r}{r_2}\right).$$

and the body force $\mathbf{f}$ as

$$\mathbf{f} = \begin{pmatrix} f_r \\ f_z \end{pmatrix} = \begin{pmatrix} 0 \\ -\frac{4\rho_{ref}u_c^2}{\rho R \Delta r} \frac{4r - \Delta r}{r} \end{pmatrix}.$$

Profiles for both the temperature and axial velocity can be found in Figure 4.

The geometry and definition of the problem are shown in Figure 5, with periodic boundary conditions on the left and right extremities of the domain. The numerical simulation is performed in a rectangular domain $\Omega: (r, z) \in [0.5, 1] \times [0, 1]$.

The convergence study is performed using a sequence of uniformly refined meshes of size $h_0 = 0.5$ for $i = 0, 1, 2, 3$, i.e. with $5 \times 10, 10 \times 20, 20 \times 40$ and $40 \times 80$ cells in the $r$ and $z$ direction respectively.

Numerical results for the Poiseuille flow between two cylinders are shown in Figure 6, in particular, the weighted $L_2$ errors for both the velocity and temperature profiles are shown for approximation orders $k = 0, 1, 2, 3, 4$, and the order of convergence of the method is found to be the optimal order $O(h^{k+1})$ for both variables.

### 4.3.2 Poiseuille flow in a tube

We now consider a slightly different configuration, where the flow is now in a tube of radius $R$ with imposed temperature $T_w$ at the wall (Figure 7).

As in the case of the flow between two cylinders, the velocity is assumed to be axial, i.e $\mathbf{u}(r, z) = u_z(r)e_z$, and that it has the following expression

$$u_z(r) = u_c \left( 1 - \frac{r^2}{R^2} \right),$$
<figure><img src="image_7.png" /><figcaption>Figure 5: Poiseuille flow between two cylinders - Problem definition.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: Convergence of the weighted $L_2$ error with the mesh size $h$ for the Poiseuille flow between two cylinders.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7: Poiseuille flow in a tube.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 8: Exact solution for the Poiseuille flow in a tube</figcaption></figure>

which is again obtained using a forcing term $f$ in the right hand side of the equations

$$ \left( \begin{array}{c} f_{\rho} \\ \boldsymbol{f} \\ f_{\rho e} \end{array} \right) = \left( \begin{array}{c} 0 \\ 0 \\ -\frac{4\mu u_c}{\rho R^2} \boldsymbol{f} \cdot \boldsymbol{u} \end{array} \right) $$

Under these conditions, an exact solution for the temperature profile $T(r)$, with $T(R) = T_w$ can be computed and is given by

$$ T(r) = T_w \left[ 1 + \frac{\text{PrMa}_{\infty}^2(\gamma - 1)}{4} \left( 1 - \frac{r^4}{R^4} \right) \right] $$

The numerical simulation is performed in the computational domain $\Omega: (r, z) \in [0, 1] \times [0, 1]$ (Figure 7b), with periodic boundary conditions imposed in the streamwise direction.

The parameters used in the numerical simulations are set as follows: we consider a flow of air, treated as a perfect gas, with heat capacity ratio $\gamma = 1.4$. The centerline velocity $u_c = 1$, the reference density $\rho_{ref} = 0.001$, the Reynolds number $Re = 1$, the wall temperature $T_w = 10$, the Mach number $Ma_{\infty} = 0.2$, and the Prandtl number $Pr = 0.7$.

The exact solution for the temperature $T(r)$, along with the imposed axial velocity profile $v(r)$, is shown in Figure 8.

Convergence of the numerical method is assessed using a sequence of quadrangular meshes with cells of size $h = \frac{h_0}{2^i}$, with $h_0 = 0.1$ and $i = 0, 1, 2, 3$, i.e. with $10 \times 10, 20 \times 20, 40 \times 40$ and $80 \times 80$ cells in the $r$ and $z$ direction respectively.

The weighted $L_2$ error is shown in Figure 9 for the axial velocity and the temperature. The order of convergence computed from these errors is the expected optimal order of convergence $O(h^{k+1})$ for both variables.

## 4.4 Flow past a sphere in the axisymmetric regime

Finally, in order to validate the implementation of the axisymmetric Navier Stokes equations, we are interested in the flow past a sphere of diameter $d$, which has been used as a benchmark test case for three dimensional codes in [13, 18, 2].
<figure><img src="image_11.png" /><figcaption>Figure 9: Convergence of the weighted $L_2$ error with the mesh size $h$ for the Poiseuille flow in a tube.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 10: Geometry, boundary conditions and mesh for the case of the axisymmetric flow past a sphere.</figcaption></figure>

In particular, despite being a rather simple configuration, the flow past a sphere can be divided into three regimes [17, 20] depending on the Reynolds number $Re_d$ (computed with respect to the diameter of the sphere $d$ and freestream quantities):

1. For small Reynolds number, namely $Re_d < 210$, the flow is expected to be in a steady axisymmetric flow regime, which is a behavior that can be observed both experimentally [30] and numerically [13, 29]. In particular, an axisymmetric recirculation zone in the wake of the sphere is observed.

2. At slightly larger Reynolds number, i.e. $210 \le Re_d \lesssim 270 - 300$, a first bifurcation happens. While the flow remains steady, the recirculation zone in the wake of the sphere is no longer axisymmetric.

3. Last, if the Reynolds is further increased to $Re_d > 270 - 300$, the flow transition into an unsteady regime and vortex shedding [28] start taking place in the wake of the sphere.

In order to validate the axisymmetric code with existing experimental and numerical data, we will focus only on the first regime at low Reynolds number.
which is the only regime for which the experimental data are axisymmetric. The code could be run on the two other regimes, but it would produce axisymmetric data which could hardly be compared to the experimental results. In particular, we will consider the flow past a sphere at Re<sub>d</sub> = 100 and Re<sub>d</sub> = 200. The test case is that of a flow of air, which is considered to be a perfect gas with specific heat ratio γ = 1.4, the gas constant R = 287 J kg<sup>-1</sup> K<sup>-1</sup>, with freestream Mach number Ma<sub>∞</sub> = 0.1, freestream pressure P<sub>∞</sub> = 100000 MPa, freestream temperature T<sub>∞</sub> = 300 K and a constant Prandtl number Pr = 0.7.

The computational domain for the axisymmetric sphere test case, along with the boundary conditions used, is shown in Figure 10a. In particular, the surface of the sphere is considered to be an adiabatic wall, and subsonic inlet and outlet conditions are used with the aforementioned freestream values.

In order to guarantee high order solutions, the surface of the sphere also needs to be meshed with curved elements (e.g. Figure 10b). Every numerical simulations in the following were run using three meshes obtained by refining (splitting) the initial mesh in Figure 10b with respectively 476, 1736, and 6608 elements.

### 4.4.1 Sphere with Re<sub>d</sub> = 100

We first consider the case of a flow past a sphere at Re<sub>d</sub> = 100, which is in the axisymmetric regime described previously. Figure 11 shows the momentum magnitude and pressure field obtained using Q³ elements with the third mesh. In particular, the recirculation zone is visible in Figure 11a, and the ratio between length of the recirculation zone ℓ and the diameter of the sphere can be computed as ℓ/d ≈ 0.875, which is around the typical value found for a sphere at Re<sub>d</sub> = 100 [18, 20].

Experimental and numerical values can be found in the literature [13, 18, 29, 23, 21] regarding, for instance, the drag coefficient C<sub>d</sub> at various Reynolds numbers. In this test case, we are interested in comparing these values to the ones obtained with the axisymmetric solver. In particular, we are interested in the convergence to expected values using high order approximations on relatively coarse meshes.

It should be emphasized that, for an axisymmetric model, integral quantities on surfaces, such as the drag coefficient C<sub>d</sub>, should be simplified in a similar way to the simplification done when deriving the numerical scheme, i.e. we have

$$C_d = \frac{2\pi}{\rho u^2 A_{\text{ref}}} \int_S ((\tau \cdot e_z) \cdot n) - P) rds,$$

where, in the case of a sphere, the reference area A<sub>ref</sub> = πd²/4.

The drag coefficient C<sub>d</sub> was computed for Re<sub>d</sub> = 100 and is shown in Table 1 and in Figure 12, and shows good agreement with the numerical and experimental data available in the literature, even on relatively coarse meshes when the approximation order is increased, except for Q⁰ computations. In particular, it can be seen that the drag coefficient converges to a value of C<sub>d</sub> ≈ 1.093, except in the case where a piecewise constant approximation Q⁰ is used. In that case, the results of the numerical simulation were not valid. For instance, the Q⁰ approximation was not able to correctly predict the recirculation zone in the wake of the sphere.

The pressure coefficient along the surface of the sphere

$$C_p(\phi) = \frac{1}{0.5\rho u_{\infty}^2} (P(\phi) - P_{\infty})$$
<figure><img src="image_13.png" /><figcaption>(a) Momentum magnitude and streamlines near the surface of the sphere, with a visible recirculation zone.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(b) Pressure isolines near the surface of the sphere. The isolines are taken as uniform values for $P \in [99500, 101000]$ with 50 samples.</figcaption></figure>

Figure 11: Numerical results for a flow past a sphere at $Re_d = 100$: Pressure and momentum fields. The results were obtained with the third mesh using $Q^3$ elements.

<table><caption>Table 1: Drag coefficient $C_d$ for the axisymmetric sphere at $Re_d = 100$. References values from Tabata et al [29]: $C_d = 1.0895$. Hartmann et al [18]: $C_d = 1.083$</caption><thead><tr><th></th><th>Q<sup>0</sup></th><th>Q<sup>1</sup></th><th>Q<sup>2</sup></th><th>Q<sup>3</sup></th><th>Q<sup>4</sup></th></tr></thead><tbody><tr><td>Mesh 1 (476 elements)</td><td>3.3983</td><td>1.4567</td><td>1.1354</td><td>1.0941</td><td>1.0928</td></tr><tr><td>Mesh 2 (1736 elements)</td><td>3.3095</td><td>1.1981</td><td>1.0991</td><td>1.0928</td><td>1.0926</td></tr><tr><td>Mesh 3 (6608 elements)</td><td>2.5603</td><td>1.1189</td><td>1.0935</td><td>1.0926</td><td>1.0926</td></tr></tbody></table>
<figure><img src="image_15.png" /><figcaption>Figure 12: Convergence of the drag coefficient</figcaption></figure>

<table><caption>Table 2: Drag coefficient $C_d$ for the axisymmetric sphere at $Re_d = 200$. References values from Tabata et al [29]: $C_d = 0.772$. Hartmann et al [18]: $C_d = 0.764$. Crivellini et al [13]: $C_d = 0.771$. Basile et al [2]: $C_d = 0.7781$.</caption><thead><tr><th></th><th>Q<sup>0</sup></th><th>Q<sup>1</sup></th><th>Q<sup>2</sup></th><th>Q<sup>3</sup></th><th>Q<sup>4</sup></th></tr></thead><tbody><tr><td>Mesh 1 (476 elements)</td><td>2.4568</td><td>1.1945</td><td>0.8247</td><td>0.7738</td><td>0.7737</td></tr><tr><td>Mesh 2 (1736 elements)</td><td>2.5977</td><td>0.8981</td><td>0.7808</td><td>0.7739</td><td>0.7739</td></tr><tr><td>Mesh 3 (6608 elements)</td><td>2.1005</td><td>0.8055</td><td>0.7749</td><td>0.7739</td><td>0.7739</td></tr></tbody></table>

where $\phi \in [0, \pi]$ is the angular position on the surface of the sphere, is shown in Figure 13 for the three meshes and shows good agreement with data from [21, 23], in particular as the order of the solution is increased.

### 4.4.2 Sphere with $Re_d = 200$

In this case, we consider a flow past a sphere with $Re_d = 200$, which is still in the axisymmetric regime. Figure 14 shows the numerical results for the momentum and pressure obtained with the third mesh and $Q^3$ elements, where the isolines and streamlines are shown for the pressure and momentum, respectively. The recirculation zone observed in Figure 14a is, as expected, longer than in the case $Re_d = 100$, and the ratio $\ell/d \approx 1.45$, which can be compared to values in the literature for $Re_d = 200$ [18, 20].

The drag coefficient $C_d$ is shown in Table 2 and in Figure 15, and it can be observed that it converges to a value $C_d = 0.7739$ when refining the mesh and approximation order, which is in agreement with values of the drag coefficient observed in the literature.

The pressure coefficient along the surface of the sphere is shown in Figure 16 for the three meshes and show good agreement with data from [21], in particular as the order of the solution is increased.
<figure><img src="image_16.png" /><figcaption>(a) First mesh</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(b) Second mesh</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(c) Third mesh</figcaption></figure>

Figure 13: Pressure coefficient $C_p(\phi)$ along the surface of the sphere at $Re_d = 100$.
<figure><img src="image_19.png" /><figcaption>(a) Momentum magnitude and streamlines near the surface of the sphere, with a visible recirculation zone.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(b) Pressure isolines near the surface of the sphere. The isolines are taken as uniform values for $P \in [99500, 101000]$ with 50 samples.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 15: Convergence of the drag coefficient</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Figure 16: Pressure coefficient $C_p(\phi)$ along the surface of the sphere at $Re_d = 200$.</figcaption></figure>
# 5 Conclusions

In this work, the development of a discontinuous Galerkin method for the simulation of flows in axisymmetric geometries was considered. In the first place, the transformation of the strong form of the Navier Stokes equations from the three dimensional form into an axisymmetric formulation using cylindrical coordinates was discussed. The formulation of the equations that was obtained resulted in additional questions about the discretization of terms which appear from geometrical manipulations on the equations, namely the non Cartesian part of the viscous flux, as well as the source term, which depends on the solution **U** and its gradient **∇U**.

Instead of starting from the equations in strong form, a discretization was derived by starting from the BR2 weak formulation of Navier-Stokes equations in three dimensions and simplifying it using three dimensional mesh induced by the rotation of a two dimensional mesh. By doing so, the aforementioned additional terms were naturally accounted for. More precisely, an expression of the numerical fluxes for both the advective and diffusive parts of the equations, as well as the discretization of the non Cartesian terms in the diffusive part were obtained. Moreover, the resulting DG method directly provides a dual consistent discretization of the geometrical source term by naturally adding the global lifting operator, which is a property required in order to ensure optimal orders of convergence in $L_2$ norm for problems with source terms depending on the gradient of the solution $\nabla U$. The axisymmetric BR2 method was then built on top of an existing code by simply adding the possibility to compute weighted integrals and adding the necessary additional terms. On the contrary to the previously proposed axisymmetric discontinuous Galerkin scheme [3], our discretization relies on a weighted lifting operator, which is consistent with the axisymmetric numerical schemes proposed in [22, 15, 8]. As our numerical scheme is obtained by averaging over the $\theta$ variable the three dimensional scheme, it inherits all the properties (stability and convergence) and the weighted $L_2$ norm, under the classical hypothesis of (9).

The accuracy of the DG discretization for axisymmetric problems was then validated on multiple and progressive test cases: scalar diffusion, advection diffusion and two types of Poiseuille flows. In these test cases, the optimal orders of convergence $O(h^{k+1})$ for approximation degrees up to $k=4$ in the weighted $L_2$ norm, which is the natural norm appearing for axisymmetric problems, were obtained. Even if our formulation and the formulation of [3] differ in the choice of the lifting operator, similar numerical results were obtained for the linear diffusion problem.

Finally, the numerical method was assessed on the problem of a laminar flow past a sphere in the axisymmetric regime, for which we were able to observe convergence of the aerodynamic coefficients (pressure, drag) to typical values found in experimental and numerical studies in the literature.

# References

[1] Douglas Norman Arnold, Franco Brezzi, Bernardo Cockburn, and L Donatella Marini. Unified analysis of discontinuous Galerkin methods for elliptic problems. *SIAM journal on numerical analysis*, 39(5):1749–1779, 2002.

[2] Francesca Basile, Jean-Baptiste Chapelier, Marta de la Llave Plata, Romain Laraufie, and Pascal Frey. Unstructured $h$-and $hp$-adaptive strategies for discontinuous Galerkin methods based on a posteriori error estimation for compressible flows. *Computers & Fluids*, 233:105245, 2022.
[3] Francesco Bassi, F Cecchi, Nicoletta Franchina, Stefano Rebay, and M Savini. High-order discontinuous Galerkin computation of axisymmetric transonic flows in safety relief valves. *Computers & Fluids*, 49(1):203-213, 2011.

[4] Francesco Bassi and Stefano Rebay. GMRES discontinuous Galerkin solution of the compressible Navier-Stokes equations. In *Discontinuous Galerkin Methods*, pages 197-208. Springer, 2000.

[5] Francesco Bassi and Stefano Rebay. Numerical evaluation of two discontinuous Galerkin methods for the compressible Navier-Stokes equations. *International journal for numerical methods in fluids*, 40(1-2):197-207, 2002.

[6] Francesco Bassi, Stefano Rebay, G Mariotti, Savini Pedinotti, and M Savini. A high-order accurate discontinuous finite element method for inviscid and viscous turbomachinery flows. In *Proceedings of the 2nd European Conference on Turbomachinery Fluid Dynamics and Thermodynamics*, pages 99-109. Antwerpen, Belgium, 1997.

[7] Yuri Bazilevs, L. Beirao da Veiga, J. Austin Cottrell, Thomas J.R. Hughes, and Giancarlo Sangalli. Isogeometric analysis: approximation, stability and error estimates for h-refined meshes. *Mathematical Models and Methods in Applied Sciences*, 16(07):1031-1090, 2006.

[8] Zakaria Belhachmi, Christine Bernardi, and Simone Deparis. Weighted Clément operator and application to the finite element discretization of the axisymmetric Stokes problem. *Numerische Mathematik*, 105:217-247, 2006.

[9] Christine Bernardi and Yvon Maday. Properties of some weighted Sobolev spaces and application to spectral approximations. *SIAM Journal on Numerical Analysis*, 1989.

[10] Anthony Bosco and Vincent Perrier. Discontinuous Galerkin methods for axisymmetric flows: extended version. Technical Report 9524, INRIA, 2023. URL: https://inria.hal.science/hal-04236220.

[11] Jean-Baptiste Chapelier, Marta De La Llave Plata, Florent Renac, and Éric Lamballais. Evaluation of a high-order discontinuous Galerkin method for the DNS of turbulent flows. *Computers & Fluids*, 95:210-226, 2014.

[12] Stephane Clain, David Rochette, and Rachid Touzani. A multislope MUSCL method on unstructured meshes applied to compressible Euler equations for axisymmetric swirling flows. *Journal of Computational Physics*, 229(13):4884-4906, 2010.

[13] Andrea Crivellini, Valerio D'Alessandro, and Francesco Bassi. Assessment of a high-order discontinuous Galerkin method for incompressible three-dimensional Navier-Stokes equations: Benchmark results for the flow past a sphere up to Re= 500. *Computers & Fluids*, 86:442-458, 2013.

[14] Simon Delmas. *Simulation d'écoulements pariétaux génériques à bas nombre de Mach pour l'amélioration du refroidissement des chambres de combustion aéronautiques*. Theses, Université de Pau et des Pays de l'Adour, December 2015. URL: https://hal.archives-ouvertes.fr/tel-01256238.

[15] Simone Deparis. *Numerical analysis of axisymmetric flows and methods for fluid-structure interaction arising in blood flow simulation*. PhD thesis, EPFL, 2004.

[16] Bruno Despres. Axisymmetric geometry. *Discontinuous Galerkin Methods: Theory, Computation and Applications*, 11:315, 2012.

[17] David Fabre, Franck Auguste, and Jacques Magnaudet. Bifurcations and symmetry breaking in the wake of axisymmetric bodies. *Physics of Fluids*, 20(5):051702, 2008.
[18] Daniel Hartmann, Matthias Meinke, and Wolfgang Schröder. A strictly conservative Cartesian cut-cell method for compressible viscous flows on adaptive grids. *Computer Methods in Applied Mechanics and Engineering*, 200(9-12):1038-1052, 2011.

[19] Ralf Hartmann. Discontinuous Galerkin methods for compressible flows: higher order accuracy, error estimation and adaptivity. *Von Karman Institute Lecture Series*, 2006.

[20] TA Johnson and VC Patel. Flow past a sphere up to a Reynolds number of 300. *Journal of Fluid Mechanics*, 378:19-70, 1999.

[21] B.P. Le Clair, A.E. Hamielec, and H.R. Pruppacher. A numerical study of the drag on a sphere at low and intermediate Reynolds numbers. *Journal of Atmospheric Sciences*, 1970.

[22] Hengguang Li. Finite element analysis for the axisymmetric Laplace operator on polygonal domains. *Journal of Computational and Applied Mathematics*, 235(17):5155-5176, 2011. URL: https://www.sciencedirect.com/science/article/pii/S0377042711002500, doi:https://doi.org/10.1016/j.cam.2011.05.003.

[23] Jacques Magnaudet, Mayela Rivero, and Jean Fabre. Accelerated flows past a rigid sphere or a spherical bubble. part 1. steady straining flow. *Journal of Fluid Mechanics*, 1995.

[24] Pierre-Henri Maire. A high-order cell-centered Lagrangian scheme for compressible fluid flows in two-dimensional cylindrical geometry. *Journal of Computational Physics*, 228(18):6882-6915, 2009.

[25] Bertrand Mercier and Genevieve Raugel. Résolution d'un problème aux limites dans un ouvert axisymétrique par éléments finis en r, z et séries de fourier en θ. *RAIRO. Analyse numérique*, 16(4):405-461, 1982.

[26] Bijan Mohammadi and Jacques Hervé Saiac. Turbulent compressible axisymmetric flows computation with the $k-ε$ model. *International Journal of Computational Fluid Dynamics*, 1(2):115-133, 1993.

[27] Todd A. Oliver and David L. Darmofal. Analysis of dual consistency for discontinuous Galerkin discretizations of source terms. *SIAM Journal on Numerical Analysis*, 47(5):3507-3525, 2009.

[28] H. Sakamoto and H. Haniu. A Study on Vortex Shedding From Spheres in a Uniform Flow. *Journal of Fluids Engineering*, 112(4):386-392, 12 1990. URL: https://doi.org/10.1115/1.2909415, arXiv:https://asmedigitalcollection.asme.org/fluidsengineering/article-pdf/112/4/386/5900198/386\_1.pdf, doi:10.1115/1.2909415.

[29] Masahisa Tabata and Kazuhiro Itakura. A precise computation of drag coefficients of a sphere. *International Journal of Computational Fluid Dynamics*, 1998.

[30] Sadatoshi Taneda. Experimental investigation of the wake behind a sphere at low Reynolds numbers. *Journal of the physical society of Japan*, 11(10):1104-1108, 1956.