Computational PDEs

Implementation of Variational Multiscale model for turbulence

The physical problem.

Starting point are the incompressible Navier-Stokes equation for a velocity \(u\) and pressure \(p\) solving \begin{align} \partial_t u - \operatorname{Re}^{-1} \Delta u + \operatorname{div}( u \otimes u ) - \nabla p & = f, \\ \operatorname{div}( u ) & = 0, \end{align} The problematic case with the incompressible Navier-Stokes equations considered here is the high Reynolds number situation which forbids the resolution of all relevant scales that may appear in a flow. In this project a Large Eddy Simulation (LES) shall be implemented, where only a certain amount of scales can be resolved, but the effect of smaller scales on the resolvable scales is modeled properly, here using a Variational Multiscale (VMS) approach.

A generic finite element discretization for the resolvable case

A generic finite element formulation transforms this into the discret formulation: Find \(u \in V_h\), \(p \in Q_h\)

\begin{align} (\partial_t u,v) + \nu a_h(u,v) + c_h(u;u,v) + b(v,p) & = f(v), \quad &&\forall v \in V_h \tag{FEMa}\\ b(u,q) & = 0, \quad &&\forall q \in Q_h \tag{FEMb} \end{align} for \(V_h\) a proper finite element space for the velocity, \(Q_h\) a proper finite element space for the pressure, \(a_h(\cdot,\cdot)\) a bilinear form for the discretization of viscosity, \(c_h(\cdot;\cdot,\cdot)\) a trilinear form discretizing convection and \(b(\cdot,\cdot)\) a bilinear form for the gradient and divergence operator discretization. In this project we want to consider the \(H(\operatorname{div})\)-conforming HDG discretization.

Variational Multiscale (VMS) modeling

Scale decomposition

In a Variational Multiscale (VMS) approach a solution \(V\) is split into three scales: large scales \(\bar{V}\), fine scales \(V'\) and unresolved scales \(V^*\). Let us denote the "big" form on the l.h.s. of (FEM) as \(K(U,V)\) where \(U=(u,p)\) and \(V=(v,q)\). The decomposition into the different scales yields the following set of equations: \begin{align} & K(\bar{U},\bar{V}) &&+ K(U',\bar{V}) &+& K(U^*,\bar{V}) &=& f(\bar{V}), \tag{VMS-I-a} \\ & K(\bar{U},V') &&+ K(U',V') &+& K(U^*,V') &=& f(V'), \tag{VMS-I-b} \\ & K(\bar{U},V^*) &&+ K(U',V^*) &+& K(U^*,V^*) &=& f(V^*). \tag{VMS-I-c} \end{align}

Getting rid of the unresolved scales (eddy-viscosity modeling)

Now, we get rid of the unresolved scale quantities \(U^*\) and \(V^*\) by the following steps:

  • by definition the unresolved scales equations (VMS-I-c) will not be solved for and will be neglected.
  • Further, the impact of the unresolved on the large scales is assumed to be zero, interaction of scales is assumed to take place through the fine scales, \(K(\bar{U},V^*)=0\).
  • Finally, the impact of the unresolved scales on the small scales is modeled through an eddy-viscosity model, i.e. a Smagorinski-type model: \(K(U^*,V') \approx \nu^* a_h(u',v')\)

We obtain: \begin{align} & K(\bar{U},\bar{V}) &&+ K(U',\bar{V}) && &=& f(\bar{V}), \tag{VMS-II-a} \\ & K(\bar{U},V') &&+ K(U',V') &+& \nu^* a_h(u',v') &=& f(V'), \tag{VMS-II-b} \\ \end{align}

Putting the scales back together (projection-based VMS)

To avoid constructing an explicit decomposition of space into two scales, we now want to put scales back together and denote now \(U = \bar{U} + U'\) and \(V = \bar{V} + V'\). Further, we introduce the large scale filter operator \(\Pi: V_h \to \bar{V}_h\) and assume the orthogonality property \(a_h(\Pi{v},v')=0\) for all \(v' \in V_h'\). We can then rewrite (VMS-II) as \begin{align} (\partial_t u,v) + (\nu+\nu^*) a_h(u,v) - \nu^* \underbrace{a_h(\Pi(u),\Pi(v))}_{=a_h(\Pi(u),v)} + c_h(u;u,v) + b(v,p) & = f(v), \quad &&\forall v \in V_h \tag{VMS-III-a}\\ b(u,q) & = 0, \quad &&\forall q \in Q_h \tag{VMS-III-b} \end{align}

A more implementation friendly formulation is finally the following: Find \(u \in V_h\), \(p \in Q_h\) and \(\bar{u} \in \bar{V}_h\) so that for all \(v \in V_h\), \(q \in Q_h\) and \(\bar{v} \in \bar{V}_h \subset V_h\) there holds

\begin{align} (\partial_t u,v) + (\nu+\nu^*) a_h(u,v) - \nu^* a_h (\bar{u},v) + c_h(u;u,v) + b(v,p) & = f(v), \quad &&\forall v \in V_h \tag{VMS-IV-a}\\ b(u,q) & = 0, \quad &&\forall q \in Q_h, \tag{VMS-IV-b}\\ \nu^* a_h(\bar{u} - u, \bar{v}) & = 0, \quad &&\forall \bar{v} \in \bar{V}_h. \tag{VMS-IV-c} \end{align} Note that here (VMS-IV-c) realizes the projection \(\Pi\) used in (VMS-III). Instead of choosing a decomposition of the space \(V_h\) explicitly we now only have to choose the large scale space \(\bar{V}_h\). Choosing \(\bar{V}_h\) large corresponds to applying the turbulence modeling on fewer scales whereas choosing \(\bar{V}_h\) as a coarse space corresponds to applying the turbulence modeling on a larger range of scales.

Another rewrite of (VMS-IV) is: Find \(u \in V_h\), \(p \in Q_h\) and \(\bar{u} \in \bar{V}_h\) such that for all \(v \in V_h\), \(q \in Q_h\) and \(\bar{v} \in \bar{V}_h\) there holds \begin{equation} (\partial_t u,v) + \nu a_h(u,v) + \nu^* a_h (u-\bar{u},v-\bar{v}) + c_h(u;u,v) + b(v,p) + b(u,q) = f(v). \tag{VMS-IV*} \end{equation}

Alternative: Projecting the stresses

Instead of applying the projection on the velocity, a similar projection can be applied on the stresses (or gradients of the velocity) yielding \begin{align} (\partial_t u,v) + (\nu+\nu^*) a_h(u,v) - 2 \nu^* \underbrace{(\Psi(\varepsilon(u)),\Psi(\varepsilon(v)))}_{=(\Psi(\varepsilon(u)),\varepsilon(v))} + c_h(u;u,v) + b(v,p) & = f(v), \quad &&\forall v \in V_h \tag{VMS-V-a}\\ b(u,q) & = 0, \quad &&\forall q \in Q_h \tag{VMS-V-b} \end{align} where \(\Psi: \nabla V_h \to \bar{S}_h\) and \(\bar{S}_h\) is a space for the large scale stresses. The more implementation friendly formulation is finally the following: Find \(u \in V_h\), \(p \in Q_h\) and \(\bar{s} \in \bar{S}_h\) so that for all \(v \in V_h\), \(q \in Q_h\) and \(\bar{r} \in \bar{S}_h\) there holds

\begin{align} (\partial_t u,v) + (\nu+\nu^*) a_h(u,v) - 2 \nu^* (\bar{s},\varepsilon(v)) + c_h(u;u,v) + b(v,p) & = f(v), \quad &&\forall v \in V_h \tag{VMS-VI-a}\\ b(u,q) & = 0, \quad &&\forall q \in Q_h, \tag{VMS-VI-b}\\ 2 \nu^* (\bar{s} - \varepsilon(u), \bar{r}) & = 0, \quad &&\forall \bar{r} \in \bar{S}_h. \tag{VMS-VI-c} \end{align} Note that here (VMS-VI-c) realizes the projection \(\Psi\) used in (VMS-V).

Most (H)DG methods \(a_h(u,v)\) can actually also be characterized through a discrete gradient operator (through lifting), cf. e.g. [LLS19b], as \(2(\varepsilon_h(u),\varepsilon_h(v))+s_h(u,v)\) for a non-negative stabilization form \(s_h(\cdot,\cdot)\), i.e. with the stresses involved in a standard \(L^2\) scalar product. A rewrite similar to (VMS-IV*) of (VMS-VI) is

Find \(u \in V_h\), \(p \in Q_h\) and \(\bar{s} \in \bar{S}_h\) such that for all \(v \in V_h\), \(q \in Q_h\) and \(\bar{r} \in \bar{S}_h\) there holds \begin{equation} (\partial_t u,v) + 2 \nu (\varepsilon_h(u),\varepsilon_h(v)) + s_h(u,v) + 2 \nu^* (\varepsilon_h(u)-\bar{s},\varepsilon_h(v)-\bar{r}) + c_h(u;u,v) + b(v,p) + b(u,q) = f(v). \tag{VMS-VI*} \end{equation} The only difference to (VMS-VI) is in the use of \(\varepsilon_h(\cdot)\) instead of \(\varepsilon(\cdot)\) in the \(\nu^\ast\)-term.

Choice for \(\nu^*\)

For \(\nu^*\) we take a Smagorinsky model with \(\nu^* = C_S \frac{h}{k+1} \Vert \varepsilon(u) \Vert_F\), where \(h\) is the local mesh size and \(k\) the (local) polynomial degree. We use \(C_S = 0.05\) as in [1].

Task

Implement and validate/evaluate the VMS method. For consistency checks you may start with 2D, but ultimately a 3D implementation and testing is necessary.

Write a solver module/class that is compatible to a standard (DNS) Navier-Stokes module/class (see NGSolve model templates). In a combination with the project ["Computational Fluid Dynamics(CFD) on HPC for (indoor air) flows"] this should allow to directly switch the HPC solver from DNS Navier-Stokes to VMS-Navier-Stokes.

References

[1] Xaver Mooslechner, Exactly Divergence-free Hybrid Discontinuous Galerkin Method for Incompressible Turbulent Flows, Diplomarbeit, TU München, 2020