A Divergence-Conforming Hybridized Discontinuous Galerkin Method for the Incompressible Reynolds Averaged Navier-Stokes Equations

TL;DR

This paper presents a novel hybridized discontinuous Galerkin method for incompressible RANS equations with turbulence modeling, achieving divergence-free velocities and efficient computation through static condensation.

Contribution

It introduces a divergence-conforming hybridized DG method that ensures point-wise divergence-free velocities and balances momentum and energy, with analysis of polynomial degree and mesh effects.

Findings

Method produces point-wise divergence-free mean velocities.

Static condensation reduces computational complexity.

Numerical results confirm effectiveness of the approach.

Abstract

We introduce a hybridized discontinuous Galerkin method for the incompressible Reynolds Averaged Navier-Stokes equations coupled with the Spalart-Allmaras one equation turbulence model. With a special choice of velocity and pressure spaces for both element and trace degrees of freedom, we arrive at a method which returns point-wise divergence-free mean velocity fields and properly balances momentum and energy. We further examine the use of different polynomial degrees and meshes to see how the order of the scalar eddy viscosity affects the convergence of the mean velocity and pressure fields, specifically for the method of manufactured solutions. As is standard with hybridized discontinuous Galerkin methods, static condensation can be employed to remove the element degrees of freedom and thus dramatically reduce the global number of degrees of freedom. Numerical results illustrate the effectiveness of the proposed methodology.

Equations

$$\Omega =[0,1{]}^2$$

$$\overline{{\bf{u}}}=\left[\aligned&-2x^{2}e^{x}(-y^{2}+y)(2y-1)(x-1)^{2}\\ -xy^{2}e^{x}(x(x+3)-2)(x-1)(y-1)^{2}\right]$$

$$\text{aligned}\overline{p}=\text{sin}(\pi x)\text{sin}(\pi y)$$

$$\text{aligned}\tilde{\nu }=\theta \text{sin}(\pi x)\text{sin}(\pi y)$$

$$\mathbf{f}=\nabla \cdot (\overline{\mathbf{u}}\otimes \overline{\mathbf{u}})+\nabla \overline{p}-\nabla \cdot (2(\nu +{\nu }_T(\tilde{\nu })){\nabla }^S\overline{\mathbf{u}})$$

$$\text{aligned}f=\overline{\mathbf{u}}\cdot \nabla \tilde{\nu }-\nabla \cdot \frac{1}{\sigma }\left[(\nu +\tilde{\nu })\nabla \tilde{\nu }\right]-C_{b1}\tilde{S}\tilde{\nu }+C_{w1}{f}_w\frac{{\tilde{\nu }}^2}{d^2}-\frac{C_{b2}}{\sigma }\nabla \tilde{\nu }\cdot \nabla \tilde{\nu }$$

Full text

1 Numerical Results

Now that we have presented our fully-discrete HDG formulation for the incompressible RANS equations coupled with the Spalart-Allmaras one equation turbulence model, we now perform verification of our method using a selection of two-dimensional example problems. We first analyze the accuracy of our formulation using the method of manufactured solutions. We then analyze the effectiveness of our formulation using three common benchmark problems: flow in a turbulent channel, flow over a backward facing step, and flow over a NACA 0012 airfoil at moderate Reynolds number and angle of attack. Each of the example problems we consider are steady in the sense that the mean velocity and pressure fields are independent of time. Nonetheless, to solve the example problems, we employed a transient analysis until a steady state flow solution was attained. All simulations were initialized using a steady state Stokes flow solution for the mean flow field and a steady state diffusion solution for the working viscosity. We also considered alternative initial conditions and found that the final steady state flow solution was independent of the initial condition in all cases.

1.1 Manufactured Solution

As a first numerical experiment, we consider a two-dimensional manufactured vortex solution for the mean flow field and a manufactured solution for the working viscosity. Namely, the flow domain for this experiment is:

[TABLE]

the mean velocity field is:

[TABLE]

the mean pressure field is:

[TABLE]

and the working viscosity is:

[TABLE]

where $\theta$ is a constant to be specified. The density $\rho$ is set to 1, and the viscosity $\nu$ is left unspecified.

The corresponding forcing in the momentum equation is:

[TABLE]

and the forcing in the Spalart-Allmaras turbulence model equation is:

[TABLE]

Homogeneous Dirichlet boundary conditions are applied along the boundary $\partial\Omega$ for both the mean velocity field and the working viscosity, and the condition $\int_{\Omega}\overline{p}d\Omega=0$ is enforced.

To assess the accuracy of our method, we have solved the manufactured solution problem for a selection of values for $\theta$ and $\nu$ as well as a series of meshes and several different polynomial degrees for the mean flow field and the working viscosity. In Fig. 4, results are displayed for $\nu=$ 1e-2, $\theta=$ 1e0, a mean velocity polynomial degree of $k$, a mean pressure polynomial degree of $k-1$, and a working viscosity polynomial degree of $k-1$ for $k=2,3,4$. From the plots, it is apparent that the mean velocity field, the mean pressure field, and the working viscosity converge at the expected rates of $k+1$, $k$, and $k$ in the $L^{2}$-norm. We observed this same behavior for other values of $\nu$ and $\theta$ as well.

To assess the impact of the polynomial degree of the working viscosity on the accuracy of the mean flow field, we also examined the impact of using a lower polynomial degree for the working viscosity than for the mean flow field. While theoretically, we expect the rate of convergence of the mean flow field to depend on the polynomial degree of both the mean flow field and the working viscosity, for certain configurations, we hypothesize that a lower polynomial degree can be employed for the working viscosity than for the mean flow field. In Figs. 8 and 12, results are displayed for $\nu=$ 1e-2, $\theta=$ 1e-2 and 1e-1 respectively, a mean velocity polynomial degree of 4, a mean pressure polynomial degree of 3, and a working viscosity polynomial degree of $q$ for $q=1,2,3$. From the plots, we observe that optimal convergence rates are attained for the mean velocity and pressure fields for $\theta=$ 1e-2, but reduced convergence rates are attained for $\theta=$ 1e-1 for higher levels of mesh refinement. These plots and further studies suggest that optimal convergence rates for the mean velocity and pressure fields are attained only pre-asymptotically, with a pre-asymptotic range whose size decreases with increasing $\theta$, unless the polynomial degree for the working viscosity is chosen to be at most one order less than the polynomial degree associated with the mean velocity. Nonetheless, for many problems of interest, including the common benchmark problems considered here, we have found that accurate mean velocity and pressure results are attainable using a low polynomial degree for the working viscosity.

1.2 Turbulent Channel Flow at $Re_{\tau}=$ 550

The next example we consider is turbulent flow in a channel at $\text{Re}_{{\tau}}$ = 550, where $\text{Re}_{{\tau}}$ is the Reynolds number based on the friction velocity and the channel half width. The domain is rectangular with dimensions 0.5 x 1.0 in the streamwise and wall normal directions, respectively. Only half the channel is modeled. For the mean flow field, periodic boundary conditions are applied in the streamwise direction, no-slip Dirichlet boundary conditions are imposed at the wall, and a symmetry condition is applied at the middle of the channel. For the working viscosity, periodic boundary conditions are applied in the streamwise direction, a homogeneous Dirichlet boundary condition is applied at the wall, and a homogeneous Neumann boundary condition is applied at the middle of the channel. 32 elements are employed in the wall-normal direction, and non-uniform grid spacing is utilized so that the boundary layer can be resolved. The polynomial degrees of the mean velocity, mean pressure, and working viscosity are taken to be 3, 2, and 1, respectively. An externally imposed pressure gradient in the streamwise direction is imposed, and the value of this gradient is selected to ensure the correct $\text{Re}_{{\tau}}$ is achieved. The kinematic viscosity is selected to be $\nu$ = 1e-4.

In Fig. 13, the mean velocity in the streamwise direction normalized by the friction velocity $u^{+}$ versus the non-dimensional distance to the wall in wall-units $y^{+}$ is displayed. Results attained using the HDG method presented in this paper are compared with those attained using a stabilized finite element implementation of the Spalart-Allmaras model within the open source CFD library PHASTA [phasta1, phasta2], as well as the DNS data of Lee and Moser [moser_direct_2015]. A grid refinement study was conducted to ensure the stabilized finite element results were fully converged. The HDG results and the stabilized finite element results are nearly indistinguishable, which is expected as they are solving the same set of partial differential equations. The HDG results also closely approximate the DNS results, which is a testament of the power of the Spalart-Allmaras model more-so than the discretization since the DNS solves a different set of partial differential equations.

1.3 Turbulent Backward Facing Step at $Re=$ 36,000

We next consider turbulent flow over a backward facing step. This problem is challenging because of the separation region that occurs downstream of the step at sufficiently large enough Reynolds numbers. The reattachment length in particular is notoriously difficult to predict. It is well-known that the Spalart-Allmaras model is not ideal for this problem, as it was created for attached wall bounded flows. However, our goal in exploring the turbulent backward facing step is not to see how well our method is able to match DNS data but rather how well it is able to match Spalart-Allmaras data in the literature and how quickly it converges with mesh refinement and polynomial degree elevation.

For all of the results reported here, a Reynolds number of $Re=$ 36,000 based on the inflow bulk velocity and step height is employed. Moreover, the kinematic viscosity is selected to be $\nu=$ 2.77777e-5, and the incoming bulk velocity is selected to be $1$. The remainder of the problem setup is as described on the NASA turbulence modeling website [nasa]. Dirichlet boundary conditions for both the mean velocity and working viscosity are applied at the inlet, with the turbulent viscosity set to be four times the kinematic viscosity. Homogeneous no-slip boundary conditions for the mean velocity and a homogeneous Dirichlet boundary condition for the working viscosity are set at the walls. Zero-traction boundary conditions for the mean velocity and a homogeneous Neumann boundary condition for the working viscosity are prescribed at the outlet.

Our first results correspond to a computational mesh consisting of 15,325 elements and mean velocity, mean pressure, and working viscosity polynomial degrees of 3, 2, and 1, respectively. The computational mesh is displayed in Fig. 14, and attained steady state solutions for the mean velocity, mean pressure, and normalized eddy viscosity are displayed in Fig. 18. In order to assess the accuracy of our results, we compare in Fig. 21 the obtained surface pressure and skin friction coefficients on the bottom surface of the domain with results publicly available on the NASA turbulence modeling website which were obtained using the NASA code CFL3D. The CFL3D results were obtained using a second-order finite volume method and a mesh of 66,049 elements, while only 15,325 elements were employed for the HDG results. Despite this discrepancy in resolution, the HDG results are visually indistinguishable from the CFL3D results. Most notably, both the HDG results and the CFL3D results predict that reattachment occurs at $x\approx 6.1$.

We next examine the impact of mesh refinement and polynomial degree elevation on the accuracy of our HDG method with a particular focus on reattachment length. In Fig. 22, the predicted reattachment length is reported for a series of refined meshes and for mean velocity, mean pressure, and working viscosity polynomial degrees of $k$, $k-1$, and 1 for $k=1,2,3$. Note that the predicted reattachment length quickly converges with mesh refinement for each of the considered polynomial degrees. Further note that the accuracy is significantly improved with polynomial degree elevation.

1.4 Turbulent Flow Over a NACA 0012 Airfoil at $\alpha$ = $10^{\circ}$

As a final numerical experiment, we analyze turbulent flow over a NACA 0012 airfoil at a Reynolds number of $Re=3,000,000$ based on the chord length and inflow velocity and an angle of attack of $\alpha=10^{\circ}$. This particular example has been explored in great depth in previous studies, and therefore there is an abundance of data to compare with for verification. We specifically setup the problem to be similar to that described on the NASA turbulence modeling website [nasa]. The domain is chosen to be a sufficiently large box about the airfoil so that the boundary conditions do not affect the flow conditions near the airfoil. The kinematic viscosity is set to $\nu=3.33333e-7$, and the mean velocity field is set to be $\overline{{\bf{u}}}=(0.9848,0.1736)$ along the left (inflow) boundary. Zero-traction boundary conditions are applied along the upper, lower, and right (outflow) boundaries, and homogeneous no-slip boundary conditions for the velocity field are applied along the airfoil surface itself. For the working viscosity, four times the kinematic viscosity is prescribed as a Dirichlet inflow condition at the left boundary, a zero Dirichlet boundary condition is applied to the surface of the airfoil, and a zero traction Neumann condition is applied at the upper, lower, and right boundaries.

To solve this flow problem, a computational mesh consisting of 61,022 elements and mean velocity, mean pressure, and working viscosity polynomial degrees of 2, 1, and 1, respectively, were employed. A conservative grid stretching ratio of approximately 1.2 was used to generate the mesh, and special care was taken to ensure that the first point off the wall has a $y+$ value of approximately 1. The computational mesh is displayed in Fig. 26, and attained steady state solutions for the mean velocity, mean pressure, and normalized eddy viscosity are displayed in Fig. 30. The attained steady state solutions match well with steady solutions appearing in the literature [burgess_high-order_2012]. To further validate our results, we compare in Fig. 33 the obtained surface pressure and skin friction coefficients on the airfoil with results publicly available on the NASA turbulence modeling website which were obtained using the NASA code CFL3D. The CFL3D results were attained using a second-order finite volume method and a mesh of 230,529 elements. The HDG and CFL3D results are visually indistinguishable, confirming the accuracy of the HDG results even though the HDG simulation only employed a quarter of the number of elements as the CFL3D simulation.