[latexpage]
[fruitful_sep]
Results
In order to verify the numerical properties of the proposed method we apply it to the solution of two well known test cases: (i) Taylor-Green vortex and (ii) inviscid shear layer roll-up. With the first test the convergence properties of the method will be assessed and with the second the time reversibility and mass, energy, enstrophy, and vorticity conservation will be verified.
1. Taylor-Green vortex
The accuracy of the proposed numerical method is compared against an analytical solution of the Navier-Stokes equations. A suitable analytical solution is the Taylor-Green vortex, given by:
\begin{equation}
\begin{cases}
u_{x} (x,y,t) = -\sin(\pi x)\cos(\pi y)\,e^{-2\pi^{2} \nu t}, \\
u_{y} (x,y,t) = \cos(\pi x)\sin(\pi y)\,e^{-2\pi^{2} \nu t}, \\
p(x,y,t) = \frac{1}{4}\,(\cos(2\pi x)+ \cos(2\pi y))\,e^{4\pi^{2} \nu t}, \\
\omega (x,y,t) = -2\pi\sin(\pi x)\sin(\pi y)\,e^{-2\pi^{2} \nu t}.
\end{cases} \label{eq:taylor_green_exact_solution}
\end{equation}
The solution is defined on the domain $\Omega = [0,2]\times[0,2]$, with periodic boundary conditions. The initial condition for both the velocity $\vec{u}(x,y,t=0)$ and vorticity $\omega(x,y,t=0)$ are given by the exact solution \ref{eq:taylor_green_exact_solution}. The kinematic viscosity is set to $\nu = 0.01$. For this study we consider the evolution of the solution from $t=0$ to $t=1$.
The first study focusses on the time convergence. For this we have used 1024 triangular elements and polynomial basis functions of degree $p=4$ and time steps equal to $\Delta t = 1, \frac{1}{2}, \frac{1}{4}, \frac{1}{8}, \frac{1}{16}$. As can be observed below, this method achieves a first order convergence rate.
Regarding the spatial convergence we have tested the convergence rate of discretizations with basis functions of different polynomial degree, $p=1,2,4$. In order not to pollute the spatial convergence rate with the temporal integration error, we have used different time steps: $\Delta t = 2.5\times 10^{-2}$ for $p=1$, $\Delta t = 1.0\times 10^{-3}$ for $p=2$ and $\Delta t = 1.0\times 10^{-4}$ for $p=4$. This method has a convergence rate of $p$-th order for basis functions of polynomial degree $p$, as can be seen below.
2. Inviscid shear layer roll-up
The second group of tests focusses on the conservation properties of the proposed numerical method. For this set of tests we consider inviscid flow $\nu=0$. We consider here the simulation of the roll-up of a shear layer. This solution is particularly challenging because during the evolution large vorticity gradients develop. Several methods are known to \emph{blow up}. We consider on the periodic domain $\Omega = [0,2\pi]\times[0,2\pi]$ the following initial conditions for velocity $\vec{u}$
\begin{equation}
u_{x}(x,y) =
\begin{cases}
\tanh\left(\frac{y-\frac{\pi}{2}}{\delta}\right), & y\leq\pi, \\
\tanh\left(\frac{\frac{3\pi}{2}-y}{\delta}\right), & y>\pi,
\end{cases}
\qquad\qquad\qquad
u_{y}(x,y) = \epsilon\sin(x),
\end{equation}
and vorticity $\omega$
\begin{equation}
\omega(x,y) =
\begin{cases}
\frac{1}{\delta}\,\text{sech}^2\left(\frac{y-\frac{\pi }{2}}{\delta }\right), & y\leq\pi, \\
-\frac{1}{\delta}\,\text{sech}^2\left(\frac{\frac{3 \pi }{2}-y}{\delta }\right), & y>\pi,
\end{cases}
\end{equation}
with $\delta = \frac{\pi}{15}$ and $\epsilon=0.05$.
The small perturbation $\epsilon$ in the $y$-component of the velocity field will trigger the roll-up of the shear layer. We show the contour lines of vorticity at $t=4$
and $t=8$
The plots on the left side correspond to a spatial discretization of 6400 triangular elements of polynomial degree $p=1$ and time step size $\Delta t = 0.1$. On the right we present the results for a spatial discretization of 6400 triangular elements of polynomial degree $p=4$ and time step size $\Delta t = 0.01$. An example mesh with 1600 elements is shown in below.
Although it is possible to observe oscillations in the vorticity plots, especially for $p=1$, none of them blows up. These oscillations necessarily appear whenever the flow generates structures at a scale smaller than the resolution of the mesh. Since energy and enstrophy are conserved, there is no possibility for the numerical simulation to dissipate the small scale information. This is a feature observed in all conserving inviscid solvers without a small scale model.
In order to assess the conservation properties of the proposed method we computed the evolution of the shear-layer problem from $t=0$ to $t=16$, using different time steps $\Delta t = 1,\frac{1}{2}, \frac{1}{4},\frac{1}{8}$. We have used a coarse grid with 6400 triangular finite elements and basis functions of polynomial degree $p=1$ for the spatial discretization. Below, on the left, we show the evolution of the kinetic energy error with respect to the initial kinetic energy, $\frac{\mathcal{K}(0)-\mathcal{K}(t)}{\mathcal{K}(0)}$. On the right, we show the evolution of the enstrophy error with respect to the initial enstrophy, $\frac{\mathcal{E}(0)-\mathcal{E}(t)}{\mathcal{E}(0)}$.
These results confirm the conservation of kinetic energy and enstrophy.
The evolution of the error of total vorticity is shown in below, on the left. On the right, we show the divergence of the velocity field.
As can be seen, this discretization conserves the total vorticity and accuratly reproduces a velocity field that is divergence free.
The final test addresses time reversibility. In order to investigate this property of the solver, we let the flow evolve from $t=0$ to $t=8$ and then reversed the time evolution and evolved again for the same time, corresponding to the evolution from $t=8$ to $t=0$. We have performed this study with 6400 triangular finite elements of polynomial degree $p=1$ and different time steps $\Delta t = 1,\frac{1}{2}, \frac{1}{4},\frac{1}{8}$. In the video below we show an animation of this simulation.
[youtube video=https://youtu.be/Su8cf5RqX5E color=red suggested=0 showinfo=0 maxw=640]