A conservative degree adaptive HDG method for transient incompressible flows

The strong form of the transient incompressible Navier-Stokes equations in an open bounded domain $\Omega\subset\mathbb{R}^{\texttt{n}_{\texttt{sd}}}$, where $\texttt{n}_{\texttt{sd}}$ is the number of spatial dimensions, is written as

where the $u$ is the velocity vector, $p$ is the pressure, $\nu$ is the kinematic viscosity, $\text{\boldmath$\nabla$\unboldmath}^{\texttt{{s}}}{\text{\boldmath$u$\unboldmath}}:=(\text{\boldmath$\nabla$\unboldmath}\text{\boldmath$u$\unboldmath}+\text{\boldmath$\nabla$\unboldmath}^{T}\text{\boldmath$u$\unboldmath})/2$ is the strain-rate tensor, $s$ is the source term, $\text{\boldmath$u$\unboldmath}_{D}$ is the imposed velocity on the Dirichlet part of the boundary, $\Gamma_{\!\!D}$, $t$ is the imposed traction on the Neumann part of the boundary, $\Gamma_{\!\!N}$, $n$ is the outward unit normal vector to the boundary, $\text{\boldmath$u$\unboldmath}_{0}$ is the initial condition and $T$ denotes the final time.

The HDG method uses a mixed formulation leading to a rewriting of the momentum equation as a first-order partial differential equation, namely

where $\text{\boldmath$L$\unboldmath}=-\sqrt{2\nu}\text{\boldmath$\nabla$\unboldmath}^{\texttt{{s}}}\text{\boldmath$u$\unboldmath}$ is the so-called mixed variable.

After discretising the domain in a set of $\texttt{n}_{\texttt{el}}$ non-overlapping elements $\Omega_{e}$, the mixed problem is written in each element and interface conditions to enforce the continuity of the solution and the continuity of the fluxes are introduced [32]. A distinctive feature of the HDG method is the introduction of the trace of the velocity, also called hybrid velocity, as an independent variable on the mesh skeleton, defined as

The resulting problem is then solved in two stages. First the element-by-element local problems define a pure Dirichlet problem

where the last equation is required to remove the indeterminacy of the pressure, $\langle\cdot,\cdot\rangle_{S}$ denotes the classical $\mathcal{L}_{2_{\!}}$ inner product for vector-valued functions on $S\subset\Gamma\cup\partial\Omega$ and $\hat{u}$ is the hybrid velocity.

Second, the global problem is given by

where the last equation is induced by the free-divergence condition in the local problems. It is worth noting that there is no need to impose the continuity of the solution in the global problem due to the unique definition of the hybrid velocity in each face of the mesh skeleton and the imposition of $\text{\boldmath$u$\unboldmath}=\text{\boldmath$\hat{u}$\unboldmath}$ in the local problems (4).

For each element, the weak formulation of local problems can be written as is as follows: find $(\text{\boldmath$L$\unboldmath}_{e},\text{\boldmath$u$\unboldmath},p)\in\left[\mathcal{H}(\operatorname{div};{{\Omega}});\SS\right]\times\left[\mathcal{H}^{1\!}(\Omega)\right]^{\texttt{n}_{\texttt{sd}}}\times\mathcal{H}^{1\!}(\Omega)$ such that

for all $(\text{\boldmath$G$\unboldmath},\text{\boldmath$w$\unboldmath},q)\in\left[\mathcal{H}(\operatorname{div};{{\Omega_{e}}});\SS\right]\times\left[\mathcal{H}^{1\!}(\Omega_{e})\right]^{\texttt{n}_{\texttt{sd}}}\times\mathcal{H}^{1\!}(\Omega_{e})$, where $\left[\mathcal{H}(\operatorname{div};{{\Omega_{e}}});\SS\right]$ denotes the space of square-integrable symmetric tensors $\SS$ of order $\texttt{n}_{\texttt{sd}}$ with square-integrable row-wise divergence.

Similarly, the weak form of the global problem reads: find $\text{\boldmath$\hat{u}$\unboldmath}\in\left[\mathcal{H}^{\frac{1}{2}}(\Gamma\cup\Gamma_{\!\!N})\right]^{\texttt{n}_{\texttt{sd}}}$ and $\text{\boldmath$\rho$\unboldmath}\in\mathbb{R}^{\texttt{n}_{\texttt{el}}}$ that satisfies

for all $\text{\boldmath$\widehat{w}$\unboldmath}\in\left[\mathcal{L}_{2_{\!}}(\Gamma\cup\Gamma_{\!\!N})\right]^{\texttt{n}_{\texttt{sd}}}$.

Following [33, 34, 32], the numerical traces appearing in the local and global problems are defined as

where $\tau^{d}$ and $\tau^{a}$ are the diffusive and convective stabilisation parameters, respectively, which in this work are defined as

where the factor 10 in the diffusive stabilisation is taken following previous work on HDG methods [25, 32] and the maximum in the convective stabilisation is taken over all the mesh nodes. Other options for the convective stabilisation, not considered here, have been recently proposed in [35].

The selected parameters ensure the satisfaction of the admissibility condition introduced in [33] to guarantee stability and well-posedness,

for some constant $\gamma$.

Introducing the numerical traces (8) into the local problems leads to the following residuals

where $\tau=\tau^{d}+\tau^{a}$. Similarly, the global problem leads to the residuals

In this work, the spatial discretisation is performed using isoparametric elements, including curved elements in the vicinity of curved boundaries. The approximation for the velocity $u$, pressure $p$, mixed variable $L$ is defined in a reference element, with polynomials of order $k\geq 1$. Similarly, the approximation of the hybrid velocity $\hat{u}$ is defined on a reference face, with polynomials of order $\hat{k}\geq 1$. The focus of this work is on degree adaptivity and, therefore, the current implementation supports an arbitrary order of approximation on each element. When two neighbouring elements use two different orders, the approximation of the hybrid velocity on the shared face between the two elements takes the maximum of the orders used on each element. This choice for the order of approximation of the hybrid velocity guarantees the optimal convergence properties of the HDG method with variable degree of approximation [36, 37].

The temporal discretisation is performed using high-order explicit first stage singly diagonal implicit Runge-Kutta (ESDIRK) integration methods. More precisely the fourth-order six-stage (ESDIRK46) scheme proposed in [38] is utilised in all the numerical examples. ESDIRK methods retain the stability properties of implicit Runge-Kutta methods and provide improved performance when compared to singly-diagonal implicit Runge-Kutta methods. In addition ESDIRK methods have been found to be more computationally efficient than other single-stage low order implicit schemes such as backward differentiation formulae (BDF) methods [39].

To strongly enforce the symmetry of the stress tensor, the present work considers the so-called Voigt notation, which has been shown [40, 25, 32] to provide the super-convergent properties described in the next section and extra efficiency when compared to formulations where the symmetry is not strongly enforced.

As usual in an HDG context [20, 41, 42, 43], hybridisation is used and in each Newton-Raphson iteration, a global problem is solved to obtain the hybrid velocity and the mean pressure, followed by the solution of multiple local problems, element-by-element, to obtain the velocity, pressure and mixed variable in each element.

For a more detailed presentation of the HDG formulation and its implementation, the reader is referred to [34, 32, 44, 45]. For a more detailed description of the Newton-Raphson linearisation strategy the reader is referred to [46, 35].

An attractive property of HDG methods is the possibility to construct a super-convergent approximation of the velocity field, also called the postprocessed velocity and denoted by $\text{\boldmath$u$\unboldmath}^{\star}$, by solving the element-by-element problem defined as

where $\tau$ is the tangential direction to the boundary.

The first equation in (13) is obtained after applying the divergence operator to the equation that defines the mixed variable and the boundary condition imposes equilibrated fluxes on the boundary of each element. The two last equations in (13) are introduced to remove the indeterminacy associated with the translational and rotational modes.

Previous work on HDG methods [47] have proved that if the velocity, pressure and mixed variable are approximated with polynomials of degree $k\geq 1$, their respective errors, measured in the $\mathcal{L}_{2_{\!}}(\Omega)$ norm, converge with order $k+1$ whereas the postprocessed velocity has an error that converges with order $k+2$, at least in diffusion dominated areas.

The possibility to build a super-convergent velocity in HDG method was first exploited in [31] to devise a cheap error indicator to drive a degree adaptive process in wave propagation problems. This strategy has also been used for incompressible Navier-Stokes flows [27], Stokes flows [48] and linear elastic problems in [49].

The main idea consists of approximating the error in the velocity field, $u$, in an element, $\Omega_{e}$, as

where the normalisation using the element measure is crucial for meshes with large variation in element size [50].

The procedure to adapt the degree of approximation aims at ensuring that the error in each element is below a user-defined tolerance $\varepsilon$ [51]. The degree is iteratively adapted as $k_{e}^{r}=k_{e}^{r-1}+\Delta k_{e}$ where $r$ denotes the degree adaptive iteration and the increment is given by

where $\lceil\cdot\rceil$ denotes the ceiling function. The base 10 in the logarithm base has been selected to minimise the number of iterations required in the degree adaptive process, but higher values can be used for a less aggressive adaptation [52, 27]. The user defined tolerance $\varepsilon$ could be selected to be a piecewise constant function with different values in different elements/regions, but for simplicity, a constant value is used in this work and detailed for each example.

For steady problems the adaptive process starts computing the solution for a given degree of approximation, commonly $k=1$ in all elements. After the solution is computed, the postprocessed velocity and the error indicator are evaluated element-by-element using (13) and (14), respectively. With this information, a new elemental degree map is defined (15). The process is repeated with the new elemental degree map until the error provided by the error indicator in each element is below the user-defined tolerance.

A solution computed with a given degree map can be used to build a better initial guess of the Newton-Raphson scheme by interpolating the solution at the new nodal distribution within each element. Let us consider that the solution in one element has been computed using a polynomial approximation of degree $r$ and the new degree to be used in the element is $s$. The solution is initially approximated as

where $\texttt{n}_{\texttt{en}}^{r}$ denotes the number of element nodes, $\mathbf{u}_{j}$ are the nodal values of the solution and $N_{j}^{r}$ are the polynomial shape functions of degree $r$ defined, on a reference element, from the set of nodes $\{\text{\boldmath$\xi$\unboldmath}^{r}\}_{i=1,\ldots,\texttt{n}_{\texttt{en}}^{r}}$. The interpolation in the new set of nodes associated to a degree $s$, $\{\text{\boldmath$\xi$\unboldmath}^{s}\}$, can be written as

where $\mathbf{u}_{j}^{s}=\text{\boldmath$u$\unboldmath}^{r}(\text{\boldmath$\xi$\unboldmath}^{s}_{j})$.

A crucial difference of a degree adaptive process for transient problems, compared to the steady case, is that the projection of the solution at time $t^{n}$ to the desired degree map is required to compute the solution at time $t^{n+1}$ and the projection is not just used as an initial guess of the Newton-Raphson scheme. Let us consider the case where the solution in one element at time $t^{n}$ is computed with a degree $r$ and the degree adaptive process changes the required degree in the element to be $s$. The projection given by (17) does not generally guarantee that the projected velocity field at time $t^{n}$ is divergence-free. More precisely, if $s\geq r$, i.e. if the adaptive process increases or maintains the degree of approximation in the element, the projection does not change the velocity field at time $t^{n}$ because the space of polynomials of degree $r$ is a subset of the space of polynomials of degree $s$. However, if $s<r$, i.e. if the adaptive process decreases the degree of approximation in the element, the projection changes the velocity field at time $t^{n}$ and the incompressibility constraint is, in general, violated.

To avoid this problem, this work proposes a new projection based on the constrained minimisation problem

The discrete version of the minimisation problem is a classical $\mathcal{L}^{2}(\Omega_{e})$ projection of the solution, whereas the constraint is imposed using a Lagrange multiplier. The resulting system of linear equations to be solved in an element where the adaptive process decreases the degree of approximation can be written as

in two dimensions, where $\mathbf{U}_{a}^{s}$ is the vector containing the nodal values of the $a$-th component of the projected free-divergence velocity field, $\lambda$ is the Lagrange multiplier,

and $u_{a}^{r}$ is the $a$-th component of the original velocity field, approximated with polynomials of degree $r$.

It is worth emphasising that the minimisation problem, i.e. the solution of the linear system (20), is only required on those elements where the adaptive process decreases the degree of approximation and the size of the linear system in two dimensions is $2\texttt{n}_{\texttt{en}}+1$ where $\texttt{n}_{\texttt{en}}$ is the number of element nodes. In addition, the problem is solved independently on each element so it can be trivially parallelised to minimise the computational overhead.

Algorithm 1 describes the degree adaptive process, including the proposed conservative projection, where $\texttt{n}_{\texttt{steps}}$ is the number of time steps, $\texttt{n}_{\texttt{adaptivity}}$ is the number of times the adaptivity is repeated each time step and $\texttt{n}_{\texttt{NR}}$ is the maximum number of iterations used in the Newton-Raphson scheme.

In the current implementation the maximum number of iterations of the Newton-Raphson scheme is five, but given the quadratic convergence only an average of three iterations are needed to reach the desired tolerance, set to $10^{-10}$ for all the residuals of the global and local problems. Given the large time steps used in the time marching process, numerical examples will be used to show that two adaptive iterations per time steps are required to obtain a converged solution with the desired error in each time step.

The first example considers the Kovasznay flow [53], which provides an analytical solution of the incompressible Navier-Stokes equations. The computational domain is a unit square, $\Omega=[0,1]^{2}$, and the analytical solution is given by

where $\lambda=\dfrac{{Re}}{2}-\sqrt{\dfrac{{Re}^{2}}{4}+4\pi^{2}}$ and $C=\dfrac{1}{8}\left[1+\exp(4\lambda)-\dfrac{1}{2\lambda}(1-\exp(4\lambda))\right]$.

A Neumann boundary condition, corresponding to the exact solution, is imposed on the bottom part of the boundary, whereas Dirichlet boundary conditions, corresponding to the exact velocity, are imposed on the rest of the boundary.

Four uniform meshes are considered, with 16, 64, 256, and 1,024 triangular elements, respectively. The first three meshes are shown in Figure 1.

Figure 2 shows the $\mathcal{L}^{2}(\Omega)$ norm of the error of the velocity, pressure, velocity gradient and postprocessed velocity as a function of the characteristic element size $h$ for a degree of approximation ranging from $k=1$ up to $k=4$.

For any degree of approximation, the expected $k+1$ convergence rate can be observed for the velocity, pressure and velocity gradient, whereas the super-convergent velocity shows the $k+2$ convergence rate. The extra accuracy of the super-convergent velocity that allows building an error indicator can be observed.

The second example, considered to verify the correct implementation of the high-order ESDIRK46 time integrator, considers the manufactured solution

where $\omega=10$ is used to define a fast variation of all flow quantities in time. The final time used in these examples is $T=0.25$ and the mesh of Figure 1(b) is used with $k=4$ to ensure that the error due to the spatial discretisation is below the error induced by the temporal discretisation.

Figure 3 shows the $\mathcal{L}^{2}(\Omega)$ norm of the error for the velocity, pressure, velocity gradient and postprocessed velocity as a function of the time step $\Delta t$.

The observed convergence rates generally align with the theoretical fourth order of accuracy for the ESDIRK46 method. The slightly lower rate observed for the pressure is associated to the so-called order reduction of ESDIRK methods [54] often observed when non-homogeneous boundary conditions are considered.

The next example considers the laminar flow, at $Re=100$, around two circular cylinders on tandem. The far field is made of a circle of diameter 100 centred at the origin, whereas the two circular cylinders have diameter 1 and are centred at $(-20,0)$ and $(10,0)$, respectively.

An unstructured mesh of 2,712 triangles is employed for this example. Curved elements are generated near the cylinder using the elastic analogy presented in [55]. Given the low Reynolds number considered, the size of the elements in the normal direction to the wall is relatively large and only the first two layers of elements around the cylinders are curved. More precisely, the size of the first element around the circular cylinders is 0.01 and the growing factor in the normal direction is 1.4. Two point sources are introduced to prescribe a mesh size of 0.2 near the cylinders, whereas a line source with size 0.75 is placed in the path of the von Karman vortex street. A detailed view of the mesh near the cylinders is shown in Figure 4.

The ESDIRK46 time marching algorithm [38] is used with a time step $\Delta t=0.2$ and the solution is advanced until the final time $T=200$.

As no analytical solution is available for this problem, a reference solution is computed by employing a uniform degree of approximation $k=6$. Further numerical experiments, not reported here for brevity, were performed to ensure that $k=6$ is the minimum degree that is required for this problem to get a converged solution. Figure 5 shows the reference pressure and magnitude of the velocity at $t=200$.

The high-order spatial approximation is crucial in this problem to accurate capture the von Karman vortex street generated by the first cylinder and its influence on the second cylinder. If a low order ($k=1$) approximation is used on the same mesh, the intensity of the vortices is not captured, as shown in Figure 6, clearly illustrating the low dissipative properties of a high-order approximation scheme.

The low order results also display a larger dispersion when compared to the high order approximation as the vortices appear in different positions.

The results shown in Figure 5 suggest that a uniform degree of approximation is not required and a degree adaptive approach is an attractive approach to increase the resolution only where is needed. The next experiments compare different degree adaptive strategies where the desired error, as detailed in Section (3.2), is taken as $\varepsilon=10^{-4}$ in all experiments, unless otherwise stated. Given the large time step used, the adaptive process is repeated twice per time step to ensure that flow features are properly captured as the solution progresses. The effect of not repeating the adaptive process is also illustrated in this example.

A standard degree adaptive approach, i.e. without the proposed correction, is first considered, where at each time step the solution in each element is projected using the desired degree of approximation map according to the error indicator provided by the HDG method. The results at $t=200$ are shown in Figure 7, including the degree used in each element.

The velocity field is in good agreement with the reference solution, with only a minor loss of intensity of the vortices behind the circular cylinders. However, the pressure field shows some important numerical artefacts when compared to the reference solution, which are related to the violation of the incompressibility constraint when projecting the velocity field from a given degree map to another degree map, in particular when the degree of approximation is decreased, as explained in Section 3.3.

To quantify the accuracy of the simulations, the lift and drag are considered the quantities of interest. Figure 8 shows the lift and drag on the first cylinder using a standard degree adaptive approach, i.e. without the proposed correction, and the result is compared to the reference solution.

The results clearly display non-physical oscillations of the drag, whereas the lift is accurately computed. Similar results for the quantities of interest for the second cylinder are shown in Figure 9.

Further experiments have been performed to confirm that the apparent more accurate results on the lift are due to the cancellation of errors and the symmetry of the lift with respect to a zero mean value. To corroborate this a mesh convergence analysis has been performed for the steady flow around a cylinder at $Re=30$, measuring the lift and drag on the upper and lower parts of the cylinders. The results show that the values of the lift and drag, measured separately on the upper and lower parts of the cylinder, converge to reference values as the mesh is refined. However, when the error of the total lift and total drag are measured, only the drag shows the expected reduction of the error as the mesh is refined, whereas the lift exhibits a very small error, even on coarse meshes, due to the addition of the upper and lower contributions, which have opposite sign.

Next, the degree adaptive procedure is enhanced by introducing the correction proposed in Section 3.3. To illustrate the benefits of the proposed approach, Figure 10 shows the degree map, pressure and magnitude of the velocity at $t=200$.

It can be observed that all the artefacts on the pressure field are not present and an excellent agreement with the reference solution is obtained.

To better quantify the accuracy of the simulation with the proposed conservative projection, Figure 11 shows the lift and drag on the first cylinder.

The results demonstrate that the proposed correction completely removes the non-physical oscillations shown in the previous simulations and provide a lift and drag which are in excellent agreement with the reference solution. The results for the second cylinder are shown in Figure 12, showing again that no oscillations are observed and a very good agreement with the reference solution is obtained.

To further illustrate the benefit of the proposed conservative projection, Table 1 reports the maximum error of the lift and drag for both cylinders.

The results clearly show the extra accuracy provided by the conservative projection. More precisely the error in the lift is more than 10 times lower using the conservative projection, whereas the error in the drag is almost 40 times lower.

To conclude this example, further numerical experiments are performed to illustrate that the conservative projection is only needed when the degree of approximation is allowed to decrease during the adaptive process. In addition, the effect of the desired error during the degree adaptive process is illustrated.

Figure 13 shows the drag on the first and second cylinders using a standard degree adaptivity where the degree of approximation is not allowed to decrease.

It can be observed that a very good agreement with the reference solution is obtained, without the oscillatory behaviour that was observed when the degree was allowed to decrease during the adaptive process. However, the main drawback of this approach is the obvious increase of computational cost because if an element reaches a high degree of approximation at one time step, the degree will be maintained at such degree for the rest of the simulation, even if there is no need to capture any features at that region for the remaining of the simulation. In this example, due to the impulsive start and the low desired error in each element $\varepsilon=10^{-4}$, all elements of the mesh require, at some instant, a degree of approximation $k=6$, so this approach leads to the same solution as the reference solution with the extra cost of computing the error indicator and projecting the solution at each time step.

If a less restrictive tolerance is used in the adaptive process, namely $\varepsilon=10^{-3}$, the quantities of interest are obtained without oscillations, as shown in Figure 14, providing evidence that the cause for the oscillations in the drag is the violation of the incompressibility condition during the projection of the solution to a lower degree.

Some discrepancies in the drag of the second cylinder are visually observed due to the use of a less restrictive tolerance.

The degree map at $t=200$ when the adaptive process is implemented without allowing the degree of approximation to be decreased and with $\varepsilon=10^{-3}$ is shown in Figure 15.

Compared to the degree map of the adaptivity process with the proposed correction, shown in Figure 10, it can be observed that the majority of elements in the wake of the two cylinders is kept to a higher degree when the adaptivity process is not allowed to lower the degree. It is also noticeable that when the degree is not allowed to decrease, a number of elements in the wake of the two cylinders end up using a degree of approximation $k=6$ whereas if the adaptivity is allowed to decrease the degree, this high degree of approximation is not required at the final time.

To quantify the reduction of degrees of freedom induced by allowing the adaptive process to decrease the degree of approximation is shown in Figure 16

The results clearly illustrate the advantage of using the proposed projection to enable the adaptive process to lower the degree during the time marching procedure. It is also worth noting that the lower the desired error, the more advantageous is to allow the degree to be lowered.

In terms of computational cost, the simulation with the proposed conservative projection is almost two times faster than the simulation with a uniform degree of approximation $k=6$. The simulation with the conservative projection is more than three times faster than the simulation not lowering the degree. The simulation not lowering the degree is actually more expensive than computing the reference solution because the majority of elements end up having the maximum degree of approximation but the cost of computing the error indicator and projecting the solution twice every time step becomes important. This shows that the reduction in degrees of freedom translates in an important reduction in computational time.

Finally, the need to repeat the adaptive process twice at each time step is also illustrated using a numerical experiment. The simulation of Figure 14 is repeated but performing the degree adaptivity only once per time step. Due to the large time step used with a high order time integrator, the computed drag shows a significant loss of accuracy, as shown in Figure 17.

The last example, inspired by [46], considers the simulation of a gust impinging on a NACA0012 aerofoil immersed in an incompressible flow at $Re=1,000$. Following [56], the gust is introduced via a localised source term. The source term in (2) is given by

where $(x_{1}^{g},0)$ denotes the centre of the rectangle of dimension $a\times b$ where the gust is generated, the wave number is given by $\alpha=\omega/v_{\infty}$ and $v_{\infty}$ the magnitude of the free-stream velocity. The constant $K$ is defined as

where $\hat{a}$ defines the region where the gust is generated, namely $\hat{a}=2\pi/a$. Finally, the functions

and

are used to guarantee a smooth transition of the flow field in the boundary of the gust region. In the current example, the parameters that define the gust are taken as $a=1$, $b=4$, $x_{1}^{g}=1.52$ and $\omega=4\pi$.

An unstructured mesh of 2,784 triangles is employed for this example. Curved elements are generated near the aerofoil using the elastic analogy presented in [55]. The size in the normal direction of the first element around the aerofoil is 0.01 and the growing factor in the normal direction is 1.2. Two point sources are introduced to prescribe a mesh size of 0.1 near the leading and trailing edges of the aerofoil, another point source is placed at the centre of the aerofoil to prescribe a size of 0.1 in the vicinity of the aerofoil, whereas a line source with size 0.4 is placed in the path of the gust. A detailed view of the mesh near the cylinders is shown in Figure 18.

Given the more complex flow dynamics of this problem, a time step $\Delta t=0.1$ and the solution is advanced using the ESDIRK46 method until a final time $T=64$. As commonly done when simulating gust around aerodynamic obstacles [56, 46] the initial condition is taken as the steady state solution of the flow around the aerofoil, in this case for $Re=1,000$. The gust is then introduced via the source term and advanced until the final time, selected so that the gust effect in the aerodynamic forces on the aerofoil disappears.

As in the previous example, a reference solution is computed by using a uniform degree of approximation $k=6$. The degree of approximation $k=6$ is selected after performing a convergence study on the fixed mesh of Figure 18. The magnitude of the velocity at some selected instants is displayed in Figure 19, showing the initial steady state solution, the perturbation of the velocity arriving and impinging on the aerofoil, the complex transient effects induced by the gust and the recovery of the steady state solution after the gust effects disappear.

The need for adaptivity in this example is even more obvious than in the previous example because the perturbation of the velocity is very localised and using a high-order approximation in the whole domain is clearly unnecessary. Next, the standard adaptive process and the adaptivity enhanced with the proposed conservative projection are considered. To remove the effect of the gust generation, when the source term that generates the gust is active, i.e. for $t\leq 10$, a constant degree of approximation $k=6$ is used in both cases. After that time the corresponding adaptive calculation is activated. This ensures that the differences in the adaptive process are not caused by a different representation of the gust. In this example, a desired error of $\varepsilon=10^{-3}$ is utilised during the adaptive process.

Figure 20 shows the lift and drag on the aerofoil using a standard degree adaptivity and the results are compared to the reference solution. As in the previous example the results show non-physical oscillations.

The oscillations are more pronounced on the drag but can also be observed on the lift in this example due to the lack of symmetry introduced by the gust. During the transient simulation, a maximum error of $2.3\times 10^{-1}$ and $3.8\times 10^{-2}$ is observed in the lift and drag respectively, clearly not providing the required accuracy for this simulation. It is worth noting that from $t=50$ to $t=51$ a constant degree of approximation, $k=6$, is used and as soon as the adaptivity is activated, a strong overshoot in the drag is observed.

When the proposed correction is introduced, an excellent agreement is again observed between the computed lift and drag and the reference solution, as shown in Figure 21.

For this example, the maximum error in the lift and drag during the whole transient process is $5.4\times 10^{-2}$ and $6.2\times 10^{-3}$, respectively, showing the extra accuracy provided by the conservative projection of the solution during the adaptive process.

To further quantify the extra accuracy provided by the proposed projection the $\mathcal{L}_{2_{\!}}([51,64])$ norm of the relative lift and drag error is computed for both adaptive approaches. Without the proposed correction the errors in lift and drag are $6.3\times 10^{-2}$ and $1.4\times 10^{-3}$ respectively, whereas when the conservative projection is used the errors in lift and drag are more than 40 times lower, namely $1.5\times 10^{-3}$ and $2.9\times 10^{-5}$.

To illustrate the ability of the degree adaptive process to accurately capture the complex flow features of this problem, lowering the degree on the elements where accuracy is no longer required, Figure 22 shows the magnitude of the velocity and the degree map at some selected instants.

Comparing the results with the reference solution of Figure 19, it can be observed that the adaptive process captures all the flow features. The degree map clearly reflects the regions where the complexity of the solution requires a higher degree of approximation to provide the desired accuracy.

In this example, the ability to lower the degree of approximation is critical to gain the benefits of a degree adaptive process, without compromising the accuracy. As the gust introduces a localised perturbation of the velocity, without lowering the degree the final degree map shows that a high order polynomial approximation is used in many areas where the flow does not show any feature. The degree map for such an approach is displayed in Figure 23.

To quantify the benefit of the proposed conservative projection, Figure 24 show the number of degrees of freedom of the global problem as a function of the non-dimensional time for the proposed approach and an adaptive process where the degree is not allowed to be decreased during the time marching process.

With the proposed projection the number of degrees of freedom at $t=64$ is 23,518 whereas for the approach not lowering the degree of approximation the number of degrees of freedom at $t=64$ reaches 45,908. The results with the conservative projection show that the most complex dynamics happen at around $t=54$, which, according to Figure 22, is precisely when the gust impinges on the aerofoil. At this point the number of degrees of freedom of the global problem reaches a maximum and then decreases because the degree of approximation can be lowered in many elements in the vicinity of the aerofoil where the transient effects are no longer relevant.

In terms of computational cost, the simulation with the proposed conservative projection is more than three times faster than the simulation with a uniform degree of approximation $k=6$. The extra performance compared to the previous example is due to the localised effect of the gust. In this example, the degree adaptive clearly offers a major advantage by introducing high order approximation only where needed.