# Hybrid RANS/DVMS modeling for static and rotating obstacles

Florian Miralles, Bastien Sauvage, Stephen F. Wornom, Bruno Koobus, Alain Dervieux

► To cite this version:

Florian Miralles, Bastien Sauvage, Stephen F. Wornom, Bruno Koobus, Alain Dervieux. Hybrid RANS/DVMS modeling for static and rotating obstacles. [Research Report] RR-9464, INRIA. 2022, pp.26. hal-03605874

HAL Id: hal-03605874

https://inria.hal.science/hal-03605874

Submitted on 11 Mar 2022

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
informatics mathematics
Inria

Modélisation hybride
RANS/DVMS autour
d'obstacles statique et
tournants

Florian Miralles, Bastien Sauvage, Stephen Wornom, Bruno Koobus, Alain Dervieux

RESEARCH
REPORT

Nº 9464

2022

Project-Team Ecuador
<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

Modélisation hybride RANS/DVMS autour
d'obstacles statique et tournants

Florian Miralles *, Bastien Sauvage †, Stephen Wornom ‡, Bruno

Koobus §, Alain Dervieux ¶

Équipe-Projet Ecuador

Rapport de recherche n° 9464 — 2022 — 26 pages

**Résumé** : Nous examinons la modélisation hybride RANS-LES pour l'étude des écoulements autour des machines tournantes comme les hélicoptères et les drones. En utilisant un maillage fin, les simulations devraient fournir des estimations précises concernant l'émission de bruit. Chacun de ces écoulements peut impliquer des régions turbulentes moyennes et hautes de Reynolds avec des tourbillons détachés et avec de fines couches limites laminaires et turbulentes. Un modèle hybride, comme le DDES, est alors obligatoire, avec éventuellement une résolution améliorée des régions LES, qui sont principalement des sillage turbulents. Il est alors intéressant d'appliquer dans ces régions un modèle LES plus sophistiqué que la partie LES de DDES. Dans notre étude, nous utilisons le modèle multi-échelle variationnel dynamique (DVMS). Dans les autres régions, une modélisation DDES ou simplement une modélisation RANS est appliquée. Dans les deux cas, une fermeture à deux équations est choisie. Après une discussion des ingrédients de la modélisation, nous présenterons une comparaison des modèles RANS, LES et hybrides pour deux séries d'écoulements. Bien que calculés par de nombreux chercheurs, les écoulements autour des cylindres restent difficiles à prévoir. La comparaison se poursuivra avec un écoulement autour d'un mélangeur en forme de croix tournant à l'intérieur d'un cylindre. Ce rapport est une version plus longue de [23]. Ce travail a été soutenu par le projet ANR NORMA dans le cadre de la subvention ANR-19-CE40-0020-01.

**Mots-clés** : modèles hybrides de turbulence, modèle variationnel multi-échelle pour LES, méthode chimera, écoulement autour d'un cylindre, machine tournante, grilles non structurées

RESEARCH CENTRE
SOPHIA ANTIPOLIS – MÉDITERRANÉE

2004 route des Lucioles - BP 93
06902 Sophia Antipolis Cedex

<footnote>* IMAG, Univ. Montpellier, CNRS, Place Eugène Bataillon, 34090 Montpellier, France. E-mail: florian.miralles@umontpellier.fr</footnote>

<footnote>† Univ. Côte d'Azur/INRIA Projet Ecuador, Sophia-Antipolis, France, E-mail: Bastien.Sauvage@inria.fr</footnote>

<footnote>‡ IMAG, Univ. Montpellier, CNRS, Montpellier, France. E-mail: stephen.wornom@inria.fr</footnote>

<footnote>§ IMAG, Univ. Montpellier, CNRS, Montpellier, France. E-mail: bruno.koobus@umontpellier.fr</footnote>

<footnote>¶ Lemma, Biot, France, and Univ. Côte d'Azur/INRIA Projet Ecuador, Sophia-Antipolis, France, E-mail: alain.dervieux@inria.fr</footnote>
Hybrid RANS/DVMS modeling for static and rotating
obstacles

**Abstract:** We examine hybrid RANS-LES modeling for the study of flows around rotating machines like helicopters and drones. Using a fine-mesh, the simulations should provide accurate estimates concerning the noise emission. Each of these flows can involve mean and high Reynolds turbulent regions with detached eddies and with thin laminar and turbulent boundary layers. A hybrid model, like DDES, is then mandatory, with possibly an improved resolution of LES regions, which are mainly turbulent wakes. It is then interesting to apply in these regions a more sophisticated LES model than the LES part of DDES. In our study, we use the dynamic variational multiscale model (DVMS). In the other regions, a DDES or simply a RANS modeling is applied. In both cases a two-equation closure is chosen. After a discussion of the modeling ingredients, we shall present a comparison of the RANS, LES, and hybrid models for two series of flows. Although computed by many researchers, flows around cylinders remain difficult to predict. The comparison will continue with a flow around a cross shaped mixing device rotating inside a cylinder. This report is a longer version of [23]. This work was supported by the ANR project NORMA under grant ANR-19-CE40-0020-01.

**Key-words:** hybrid turbulence models, variational multiscale model for LES, chimera method, cylinder flow, rotating machine, unstructured grids
Table des matières

<table><tr><td>1</td><td>Introduction</td><td>5</td></tr><tr><td>2</td><td>Turbulence modeling</td><td>5</td></tr><tr><td>2.1</td><td>RANS</td><td>5</td></tr><tr><td>2.2</td><td>DDES model</td><td>7</td></tr><tr><td>2.3</td><td>High Reynolds computation</td><td>8</td></tr><tr><td>2.4</td><td>LES-like model : DVMS</td><td>8</td></tr><tr><td>2.4.1</td><td>SGS model</td><td>9</td></tr><tr><td>2.4.2</td><td>DVMS</td><td>10</td></tr><tr><td>2.5</td><td>Hybrid RANS/DVMS</td><td>12</td></tr><tr><td>2.6</td><td>Hybrid DDES/DVMS</td><td>13</td></tr><tr><td>3</td><td>Numerical discretization</td><td>13</td></tr><tr><td>3.1</td><td>Numerical scheme</td><td>13</td></tr><tr><td>3.2</td><td>Mesh adaptation for rotating machines</td><td>15</td></tr><tr><td>4</td><td>Applications</td><td>16</td></tr><tr><td>4.1</td><td>Flow past a cylinder</td><td>16</td></tr><tr><td>4.2</td><td>Flow around a rotating cross</td><td>23</td></tr><tr><td>5</td><td>Conclusion</td><td>23</td></tr></table>
NOMENCLATURE

<table><tr><td>W</td><td>[-]</td><td>flow variables</td></tr><tr><td>ρ</td><td>[kg/m³]</td><td>density</td></tr><tr><td>u</td><td>[m/s]</td><td>velocity vector</td></tr><tr><td>E</td><td>[J/m³]</td><td>total energy per unit volume</td></tr><tr><td>k</td><td>[m²/s²]</td><td>turbulence kinetic enery</td></tr><tr><td>ε</td><td>[m²/s³]</td><td>dissipation rate of k</td></tr><tr><td>W<sub>h</sub></td><td>[-]</td><td>discrete flow variables</td></tr><tr><td>W'<sub>h</sub></td><td>[-]</td><td>small resolved scales of W<sub>h</sub></td></tr><tr><td>⟨⟩</td><td>[-]</td><td>Reynolds average</td></tr><tr><td>F</td><td>[-]</td><td>convective and diffusive fluxes</td></tr><tr><td>τ<sup>RANS</sup></td><td>[-]</td><td>RANS closure term</td></tr><tr><td>τ<sup>LES</sup></td><td>[-]</td><td>LES closure term</td></tr><tr><td>τ<sup>DDES</sup></td><td>[-]</td><td>DDES closure term</td></tr><tr><td>ν</td><td>[m²/s]</td><td>kinematic viscosity of the fluid</td></tr><tr><td>ν<sub>t</sub></td><td>[m²/s]</td><td>turbulent kinematic viscosity</td></tr><tr><td>C&#773;<sub>d</sub></td><td>[-]</td><td>mean drag coefficient</td></tr><tr><td>C'<sub>l</sub></td><td>[-]</td><td>r.m.s. of the lift coefficient</td></tr><tr><td>C<sub>pb</sub></td><td>[-]</td><td>mean base pressure coefficient</td></tr><tr><td>θ&#773;</td><td>[deg.]</td><td>mean separation angle</td></tr><tr><td>Δ</td><td>[m]</td><td>LES filter width</td></tr><tr><td>l<sub>RANS</sub></td><td>[m]</td><td>RANS turbulence length scale</td></tr><tr><td>μ<sub>SGS</sub></td><td>[kg/m/s]</td><td>LES eddy viscosity</td></tr><tr><td>μ<sub>RANS</sub></td><td>[kg/m/s]</td><td>RANS eddy viscosity</td></tr><tr><td>τ&#771;</td><td>[-]</td><td>Reynolds stress tensor</td></tr></table>

Sets

<table><tr><td>Ω<sub>f</sub></td><td>Fluid domain</td></tr><tr><td>Ω<sub>h</sub></td><td>Collection of Tetrahedron</td></tr><tr><td>Ξ<sub>k</sub></td><td>Index set which are contained in macro cell C<sub>k</sub></td></tr><tr><td>C<sub>k</sub></td><td>Agglomeration cell</td></tr><tr><td>V<sub>2h</sub></td><td>the space of large scales</td></tr><tr><td>V'<sub>h</sub></td><td>the space of small scales</td></tr></table>
# 1 Introduction

In this paper, hybrid models are evaluated for the simulation of massively separated flows around fixed and moving geometries, with the objective of applying them to aeroacoustic problems involving complex industrial flows at high Reynolds numbers. For this purpose, the turbulence models need to behave preferentially like large eddy simulation models (LES) in order to introduce little dissipation and thus better capture the small scales and the fluctuations of the unsteady flows considered. However, these flows involve very thin boundary layers, which current computers and softwares cannot compute in LES mode, only in RANS mode. To accommodate both needs, LES and RANS, the hybrid turbulence approach, combining LES and RANS in a somewhat zonal manner, is considered by many research teams.

In this work, beside the classical Delayed Detached Eddy Simulation model (DDES, [32]), we study an approach which combines zonally this last model and the dynamic variationnal multiscale model (DVMS, [19]). In this hybrid DDES/DVMS approach [24], the DVMS model is preferentially activated in the wake regions where this latter model introduces less dissipation than the LES component of DDES. A hybrid RANS/DVMS model [16], where the RANS component is Goldberg's $k-\epsilon$ model [10], is also applied in this study. A smooth blending function, which is based on the value of a blending parameter, is employed in these hybrid strategies, for switching from RANS to DVMS or DDES to DVMS. In [16], these models are applied to flows around a cylinder at moderate Reynolds number, starting from Reynolds number 20000.

In many industrial cases, and particularly cases with rotating geometries like propellers, the Reynolds number is much higher.

In the present paper, these hybrid models are applied to the flow around a circular cylinder in the supercritical regime, namely at Reynolds numbers 1 million and 2 million. At Reynolds number 1 million, the supercritical regime shows turbulent flow separation on both the upper and lower surface of the cylinder with a laminar-turbulent transition in the boundary layers located between the stagnation point and the separation one. At Reynolds number 2 million, the high transition regime shows a fully turbulent boundary layer on one side of the cylinder, and a partially laminar, partially turbulent boundary layer on the other side. This benchmark, which contains many characteristics encountered in industrial flows, is challenging due to the complex physics of the flow and the considered high Reynolds numbers, and remains a difficult calculation to perform. The second application concerns the flow at Reynolds number 1.8 million around a cross shaped mixing device rotating inside a cylinder.

This work is part of a cooperation program with Keldysh Institute of Applied Mathematics of Moscow focussing on the "Efficient simulation of noise of rotating machines", see e.g. [7, 3, 6].

The remainder of this paper is organized as follows. Section 2 presents the hybrid turbulence models used in this work. Section 3 describes the numerical scheme. In Section 4, applications are presented. The obtained results are analyzed and compared with those obtained in other numerical studies and with experimental data. Finally, concluding remarks are drawn in Section 4.

# 2 Turbulence modeling

## 2.1 RANS

First, we want to specify that RANS stands for unsteady RANS throughout this document. When comparing with other published results, we use the notation that appear in the published
articles notable URANS which is shorthand for Unsteady RANS, which the same as RANS in this document. Steady RANS will be denoted as SRANS and Steady LES as SLES.

In this work, and as far as the closure of the RANS equations is concerned, we use the low Reynolds $k - \epsilon$ model proposed in Goldberg et al.[10]. This model was designed to improve the predictions of the standard $k - \epsilon$ one for adverse pressure gradient flows, including separated flows. The equations of this model are recalled in this document, in dimensionless frame :

$$ \mu_t = \bar{\rho} c_\mu f_\mu(k, \epsilon) \frac{k^2}{\epsilon}, \qquad (1) $$

where the damping function is defined as follow :

$$ f_{\mu}(k, \epsilon) = \frac{1 - e^{a_{\mu} R_t(k, \epsilon)}}{1 - e^{-\sqrt{R_t(k, \epsilon)}}} \max \left( 1, \frac{c_{\tau}}{\sqrt{R_t(k, \epsilon)}} \right), \quad R_t(k, \epsilon) = Re \frac{k^2}{\epsilon}. \qquad (2) $$

function $R_\mu$ represent the turbulent Reynolds number. Transports equation on $k$ and $\epsilon$ can be written as :

$$ \frac{\partial \bar{\rho} k}{\partial t} + \nabla \cdot (\bar{\rho} \tilde{\mathbf{u}} k) = \tilde{\tau} : \nabla \tilde{\mathbf{u}} + \nabla \cdot \left[ \left( \frac{1}{Re} + \frac{\mu_t}{\sigma_k} \right) \nabla k \right] - \bar{\rho} \epsilon, \qquad (3) $$

$$ \frac{\partial \bar{\rho} \epsilon}{\partial t} + \nabla \cdot (\bar{\rho} \tilde{\mathbf{u}} \epsilon) = \nabla \cdot \left[ \left( \frac{1}{Re} + \frac{\mu_t}{\sigma_\epsilon} \right) \nabla \epsilon \right] + \left( c_\epsilon^{(1)} \frac{\epsilon}{k} \tilde{\tau} : \nabla \tilde{\mathbf{u}} - c_\epsilon^{(2)} \bar{\rho} \frac{\epsilon^2}{k} + C^{(2)}(k, \epsilon) \right) C^{(1)}(k, \epsilon). \qquad (4) $$

Since to control production and dissipation of dissipation rate, functions $C^{(1)}$ and $C^{(2)}$ are introduced. These functions are defined such that :

$$ \begin{aligned} C^{(1)}(k, \epsilon) &= \frac{k}{\epsilon} \max \left( 1, \frac{c_{\tau}}{\sqrt{R_t}} \right), \\ C^{(2)}(k, \epsilon) &= \frac{3}{10} \bar{\rho} \max \left( \sqrt{k}, \left( \frac{1}{Re} \epsilon \right)^{1/4} \right) \sqrt{\epsilon C^{(1)}} \max \left( \nabla k \cdot \nabla \left( \frac{k}{\epsilon} \right), 0 \right). \end{aligned} $$

The total system, including the fluids variables and turbulent variables, can be written as :

$$ \left\{ \begin{array}{l} \frac{\partial W}{\partial t} + \nabla \cdot \underline{F}(W) = \tau^{\text{RANS}}(W), \\ + \text{Boundary conditions}, \\ Q(0, \mathbf{x}) = Q_0(\mathbf{x}), \quad \forall \mathbf{x} \in \Omega_f \end{array} \right. \qquad (5) $$

Where $\underline{\tau}^{\text{RANS}}$ means the RANS closure terms, defined such that :

$$ \tau^{\text{RANS}}(W) = \left( \underbrace{\rho \mathbf{u}}_{0}, \underbrace{\rho \mathbf{u}}_{0}, \underbrace{\rho E}_{0}, \overbrace{\rho k}^{\rho \epsilon}, \overbrace{\left( c_{\epsilon}^{(1)} \frac{\epsilon}{k} \tilde{\tau} : \nabla \tilde{\mathbf{u}} - c_{\epsilon}^{(2)} \bar{\rho} \frac{\epsilon^2}{k} + C^{(2)}(k, \epsilon) \right) C^{(1)}(k, \epsilon)}^{\rho \epsilon} \right) \qquad (6) $$

In practice, for compressible flow, we use the constants summarize in the table below :

<table><thead><tr><th>c<sub>ε</sub><sup>(1)</sup></th><th>c<sub>ε</sub><sup>(2)</sup></th><th>σ<sub>k</sub></th><th>σ<sub>ε</sub></th></tr></thead><tbody><tr><td>1.44</td><td>1.92</td><td>1.0</td><td>1.3</td></tr></tbody></table>
## 2.2 DDES model

The DDES approach used in this work is based on the Goldberg's $k - \epsilon$ model [10]. In other words, this RANS model is introduced in the DDES formulation [32] by replacing the $D_k^{RANS} = \rho \epsilon$ dissipation term in the $k$ equation by $D_k^{DDES} = \rho l_{DDES}^{k^{3/2}}$ where

$$l_{DDES} = \frac{k^{3/2}}{\epsilon} - f_d \times \max(0, \frac{k^{3/2}}{\epsilon} - C_{DDES}\Delta)$$
with $f_d = 1 - \tanh((8r_d)^3)$
and $r_d = \frac{\nu_t + \nu}{\max(\sqrt{u_{i,j}u_{i,j}}, 10^{-10})K^2d_w^2}$.

$K$ denotes the von Karman constant ($K = 0.41$), $d_w$ the wall-normal distance, $u_{i,j}$ the $x_j$-derivative of the $i$th-component of the velocity $u$, and the model constant $C_{DDES}$ is set to the standard value 0.65. We define here the DDES $k - \epsilon$ Goldberg source term :

$$\underline{\tau}^{DDES}(\underline{W}) = \left( \underbrace{0}_{\rho}, \underbrace{\rho\mathbf{u}}_{\rho\mathbf{u}}, \underbrace{0}_{\rho E}, \underbrace{\rho^k}_{\tilde{\tau}: \nabla\tilde{\mathbf{u}} - D_k^{DDES}}, \overbrace{\left( c_{\epsilon}^{(1)}\frac{\epsilon}{k}\tilde{\tau}: \nabla\tilde{\mathbf{u}} - c_{\epsilon}^{(2)}\bar{\rho}\frac{\epsilon^2}{k} + C^{(2)}(k, \epsilon) \right) C^{(1)}(k, \epsilon)}^{\rho\epsilon} \right) \quad (7)$$

Symbol $\Delta$ holds for the filter width, in the type of LES which we use, the filter width is fixed to the local mesh size. The definition of the local mesh size is not an easy question when (highly) anisotropic meshes are used. Two main options are (1) the third root of the grid element volume and (2) the largest local edge. However, our experience is that, when combined with devices like VMS and Dynamic-LES, the choice in defining the local mesh size is not so sensitive. Our filter width is defined as the third root of volume of the grid element T, then :

$$\Delta = \Delta_T = \left( \int_T d\mathbf{x} \right)^{1/3} . \qquad (8)$$
## 2.3 High Reynolds computation

In order to catch the turbulent boundary layer, we use a wall law modeling, designed by Reichardt [27]. By this way, the no slip boundary condition on the wall become:

$$U^+ = \frac{1}{\kappa} \ln (1 + \kappa y^+) + 7.8 \left[ 1 - \exp \left( \frac{-y^+}{11} \right) - \frac{-y^+}{11} \exp \left( \frac{-y^+}{3} \right) \right] \quad (9)$$

Also, for turbulent variables, we use the work of Jaeger [17], with the following condition on the wall:

$$k = \frac{u_f^2}{\sqrt{C_\mu}} \left(\frac{y^+}{10}\right)^2, \\ \varepsilon = Re \frac{u_f^4}{10\kappa} \left( \left(\frac{y^+}{10}\right)^2 + \frac{0.2}{\sqrt{C_\mu}} \left[1 - \left(\frac{y^+}{10}\right)^2\right] \right).$$

The Reichardt wall law give us a smooth transition from linear to logarithmic boundary layer.

## 2.4 LES-like model : DVMS

In this section, we present the DVMS model which is preferred to the classical LES approach in our hybrid strategies because of some specific properties, as will be specified hereafter, that allow this model to have a better behavior in some regions of the flow, such as shear layers and wakes.

The Variational Multiscale (VMS) model for the large eddy simulation of turbulent flows has been introduced in [15] in combination with spectral methods. In [19], an extension to unstructured finite volumes is defined. This method is adapted in the present work. Let us explain this VMS approach in a simplified context. Assume the mesh is made of two embedded meshes. On the fine mesh we have a $P^1$-continuous finite-element approximation space $V_h$ with the usual basis functions $\Phi_i$ vanishing on all vertices but vertex $i$. Let be $V_{2h}$ its embedded coarse subspace $V_{2h}$. Let $V'_h$ be the complementary space : $V_h = V_{2h} \oplus V'_h$. The space of small scales $V'_h$ is spanned by only the fine basis functions $\Phi'_i$ related to vertices which are not vertices of $V_{2h}$. Let us denote the compressible Navier-Stokes equations by : $\frac{\partial W}{\partial t} + \nabla \cdot F(W) = 0$ where $W = (\rho, \rho u, E)$.

The VMS discretization writes for $W_h = \sum W_i \Phi_i$:

$$\left\{ \begin{array}{l} \left( \frac{\partial W_h}{\partial t}, \Phi_i \right) + \left( \nabla \cdot F(W_h), \Phi_i \right) = - \left( \bar{\tau}^{LES}(W_h'), \Phi_i' \right) \quad \forall \mathbf{a}_i \in \Omega_h, \\ Q_h(0, \mathbf{a}_i) = Q_0(\mathbf{a}_i), \quad \forall \mathbf{a}_i \in \Omega_h \end{array} \right. \quad (10)$$

Here the convective and viscous flux do not contains any closure equations. Moreover, the closure terms is defined as the following:

$$(\bar{\tau}^{LES}(W_h'), \phi_i') := (0, \mathbf{M}_S(W_h', \Phi_i'), \mathbf{M}_H(W_h', \Phi_i'))^T \quad (11)$$

and

$$\mathbf{M}_S(W_h', \Phi_i') := \sum_{T \in \Omega_h} \int_T \mu_i^{sgs} P \nabla \Phi_i' d\mathbf{x} \quad (12)$$

Inria
here $P = 2S - \frac{1}{3}Tr(S)Id$. The basis function are written as $\Phi'_k = \Phi_k - \bar{\Phi}_k$ where :

$$ \bar{\Phi}_k = \frac{Vol(C_k)}{\sum_{j \in \mathcal{J}_k} Vol(C_j)} \sum_{j \in \mathcal{J}_k} \Phi_j \qquad (13) $$

$\mathcal{J}_k$ means the index set of cells which are contained in macro-cell $\mathcal{C}_k$ obtained by agglomeration process, see also [20]. For the subgrid terms in relation with heat transfer, we write :

$$ M_H(\underline{W}_h, \Phi'_i) = \sum_{T \in \Omega_h} \int_T \frac{C_p \mu_t^{sgs}}{Pr_t} \nabla T' \cdot \nabla \Phi'_i d\mathbf{x} \qquad (14) $$

For a test function related to a vertex of $V_{2h}$, the RHS vanishes, which limits the action of the LES term to small scales. In practice, embedding two unstructured meshes $V_h$ and $V_{2h}$ is a constraint that we want to avoid. The coarse level is then built from the agglomeration of vertices/cells as sketched in Fig. 1. It remains to define the model term $\underline{\tau}^{LES}(\underline{W}'_h)$. This term

<figure><img src="image_4.png" /><figcaption>FIGURE 1 – Building the VMS coarse level</figcaption></figure>

represents the subgrid-scale (SGS) stress term, acting only on small scales $W'_h$, and computed from the small scale component of the flow field by applying either a Smagorinsky [31] or a WALE SGS model [26], the constants of these models being evaluated by the Germano-Lilly dynamic procedure [8, 21]. The resulting model is denoted DVMS in this paper. It has been checked [25] that combining VMS and dynamic procedure effectively brings improved predictions.

A key property of the VMS formulation is that the modeling of the dissipative effects of the unresolved structures is only applied on the small resolved scales, as sketched in Fig. 2. This property is not satisfied by LES models which damp also the other resolved scales. Important consequences are that a VMS model introduces less dissipation than its LES counterpart (based on the same SGS model) and that the backscatter transfer of energy from smallest scales to large scales is not damped by the model. These VMS models then generally allow a better behavior near walls, in shear layers and in the presence of large coherent structures.

### 2.4.1 SGS model

To obtain the SGS viscosity needed to close the above VMS formulation, we have to define the SGS viscosity coefficient $\mu_t^{sgs}$. A first option is the widely used Smagorinsky model [31] is first considered. In the adopted VMS formulation this writes :

$$ \mu_t^{sgs} = \rho(C_s \Delta)^2 |\tilde{S}|, \quad \text{with } |\tilde{S}| = \sqrt{\tilde{S} : \tilde{S}}, \qquad (15) $$

where $C_s$ is the Smagorinsky coefficient. A typical value for the Smagorinsky coefficient for shear flows is $C_s = 0.1$, which is used herein.
A second SGS model is the Wall-Adapting Local Eddy -Viscosity (WALE) SGS model proposed by Nicoud and Ducros [26]. The eddy-viscosity term in the VMS formulation is defined as follows :

$$ \mu_t^{sgs} = \rho(Cw\Delta)^2 \frac{(\tau^W : \tau^W)^{3/2}}{(\check{S} : \check{S})^{5/2} + (\tau^W : \tau^W)^{5/2}}. \quad (16) $$

with the tensor invariant :

$$ \tau_{i,j}^{W} = \frac{1}{2} (\check{g}_{i,j}^2 + \check{g}_{j,i}^2) - \frac{1}{2} \delta_{i,j} \check{g}_{k,k}^2, \quad \check{g}_{i,j}^2 := \sum_k \frac{\partial \tilde{u}_i}{\partial \mathbf{x}_k} \frac{\partial \tilde{u}_i}{\partial \mathbf{x}_k}. \quad (17) $$

In this model Cw is set as 0.5, as indicated [26].

## 2.4.2 DVMS

Moreover, in this work, the dynamic procedure, which provides a tuning of the SGS dissipation in space and time, is combined with the VMS approach, which limits its effects to the smallest resolved scales, so that the resulting DVMS model enjoys synergistic effects.

<figure><img src="image_5.png" /><figcaption>FIGURE 2 – VMS principle.</figcaption></figure>

In their original formulations, Cs and Cw appearing in the expression of the viscosity of the Smagorinsky and WALE SGS model (Eqs. 15 and 16 respectively) are set to a constant over the entire flow field and in time. In the dynamic model [8], this constant is replaced by a dimensionless parameter C(x, t) that is allowed to be a function of space and time. The dynamic approach also provides a systematic way for adjusting this parameter in space and time by using information from the resolved scales. After the introduction of the grid filter, denoted by overline and tilde, tilde being Favre averaging, $\tilde{f} = \frac{\rho \tilde{f}}{\tilde{\rho}}$, a second filter is considered, having a larger width than the grid one, which is called the test-filter and denoted by a hat. The test-filter is applied to the grid-filtered Navier-Stokes equations, and then, the subtest-scale stress is defined as follows :

$$ \mathcal{T}^{\text{test}} = \widehat{\rho \mathbf{u} \otimes \mathbf{u}} - \frac{\widehat{\rho \mathbf{u}} \otimes \widehat{\rho \mathbf{u}}}{\widehat{\rho}}, \quad \text{and} \quad \tau^{\text{sgs}} = \overline{\rho \mathbf{u} \otimes \mathbf{u}} - \overline{\rho \tilde{\mathbf{u}} \otimes \tilde{\mathbf{u}}} \quad (18) $$

Then we consider the resolved turbulent tensor :

$$ \mathcal{L} = \widehat{\rho \tilde{\mathbf{u}} \otimes \tilde{\mathbf{u}}} - \frac{\widehat{\rho \mathbf{u}} \otimes \widehat{\rho \mathbf{u}}}{\widehat{\rho}} \quad (19) $$

This tensor can be interpreted as an intermediate grid filter tensor. We related these quantity by the following equation :

$$ \mathcal{L} = \mathcal{T} - \hat{\tau}, \quad (20) $$
By Smagorinsky hypothesis [31], the deviatoric part of subgrid scale tensor is proportional to stress tensor $\tilde{P} = 2\tilde{S} - \frac{2}{3}Tr(S)Id$:

$$ \text{dev } \tau = \tau - \frac{1}{3} \text{Tr}(\tau) \text{Id} = -\bar{\rho}(C_s \Delta)^2 g(\tilde{\mathbf{u}}) \tilde{P}. \quad (21) $$

$$ \text{dev } \mathcal{T} = \mathcal{T} - \frac{1}{3} \text{Tr}(\mathcal{T}) \text{Id} = -\hat{\bar{\rho}}(C_s \tilde{\Delta})^2 \hat{g}(\tilde{\mathbf{u}}) \hat{\tilde{P}} \quad (22) $$

and where $g(\tilde{\mathbf{u}})$ denotes the contribution to the SGS viscosity depending on the gradient velocity that appears in 15, explicitly $g(\tilde{\mathbf{u}}) = |\tilde{S}|$ for the Smagorinsky model, and in 16, explicitly $g(\tilde{\mathbf{u}}) = \frac{(\tau^{W;\tau W})^{3/2}}{(\tilde{S}:\tilde{S})^{5/2}+(\tau^{W;\tau W})^{5/2}}$, for the WALE model.

Then, from the resolved field, we can explicitly compute $\mathcal{L}$. And, moreover $\mathcal{L} = \text{dev}\mathcal{L}$, we have:

$$ \text{dev}\mathcal{L} = \mathcal{L} - \frac{1}{3}\text{Tr}(\mathcal{L})\text{Id}, \quad (23) $$

$$ = \mathcal{T} - \hat{\tau} - \frac{1}{3}\text{Tr}(\mathcal{T} - \hat{\tau})\text{Id}, \quad (24) $$

$$ = \mathcal{T} - \frac{1}{3}\text{Tr}(\mathcal{T})\text{Id} - \left(\hat{\tau} - \frac{1}{3}\text{Tr}(\hat{\tau})\text{Id}\right), \quad (25) $$

$$ = -(C\tilde{\Delta})^2 \hat{\bar{\rho}}g(\tilde{\mathbf{u}})\hat{\tilde{P}} + (C\Delta)^2 \hat{\bar{\rho}}g(\tilde{\mathbf{u}})\hat{\tilde{P}} \quad (26) $$

$$ = (C\Delta)^2 \underbrace{\left( \hat{\bar{\rho}}g(\tilde{\mathbf{u}})\hat{\tilde{P}} - \left( \frac{\hat{\Delta}}{\Delta} \right)^2 \hat{\bar{\rho}}g(\tilde{\mathbf{u}})\hat{\tilde{P}} \right)}_{\mathcal{B}}. \quad (27) $$

Therefore, the problem becomes to find $C\Delta$ solution of equation below:

$$ \mathcal{L} - (C\Delta)^2\mathcal{B} = 0 \Leftrightarrow \min \|\mathcal{L} - (C\Delta)^2\mathcal{B}\|^2. \quad (28) $$

By convexity and a first optimality order, the least square problem provide a uniqueness solution:

$$ (C\Delta)^2 = \frac{\mathcal{B} : \mathcal{L}}{\mathcal{B} : \mathcal{B}}. \quad (29) $$

In practice this filter is chosen such that $\hat{\Delta} > 1$. In our case $\Delta$ is defined by:

$$ \Delta = \left( \frac{1}{N_{a_i}} \int_{T_{a_i}} dx \right)^{1/3}, \\ \hat{\Delta} = \left( \int_{T_{a_i}} dx \right)^{1/3}. $$

Thus $(\hat{\Delta}/\Delta)^2 = (N_{a_i})^{2/3}$ greater than one.

Note that all quantities in the right-hand side of Eq. 29 are known from the LES computation. Note also that we preferred to dynamically compute $(C\Delta)^2$, instead of $C$ as done in the original dynamic procedure, in order to partially overcome difficulties in the definition of the filter width for inhomogeneous and unstructured grids. Finally, as done also in [14], the classical dynamic
procedure previously briefly outlined, which involves all the resolved scales, is used herein. Once $(C\Delta)^2$ is dynamically computed, it is injected in Eq. 15 or 16 to obtain the SGS viscosity used in the VMS approach.

A possible drawback of the dynamic procedure based on the Germano-identity [8] when applied to a SGS model already having a correct near-wall behavior, as the WALE one, is the introduction of a sensitivity to the additional filtering procedure. A simple way to avoid this inconvenient is to have a sensor able to detect the presence of the wall, without a priori knowledge of the geometry, so that the dynamic SGS model adapts to the classical constant of the model, which is equal to 0.5 in the near wall region for the WALE model, and compute the constant dynamically otherwise. We adopt the sensor proposed in [36], having the following expression :

$$f_{SVS} = \frac{(\tau^W : \tau^W)^{3/2}}{(\tau^W : \tau^W)^{3/2} + (\tilde{S} : \tilde{S})^3} \quad (30)$$

This parameter has the properties to behave like $y^+$ near a solid wall, to be equal to 0 for pure shear flows and to 1 for pure rotating flows. It should be noticed that the implementation of the dynamic SGS models in our software has been optimized so that the additional cost of the resulting dynamic LES and VMS models, in the case of an implicit time-marching scheme, which is our default option, is less than 1% compared to their non-dynamic counterparts.

## 2.5 Hybrid RANS/DVMS

The central idea of the hybrid RANS/DVMS model [24] applied in this study is to combine the mean flow field obtained by the RANS component with the application of the DVMS model wherever the grid resolution is adequate. In this hybrid model, the $k-\epsilon$ model of Goldberg (subsection 2.1) is used as the RANS component. First, let us write the semi-discretization of the RANS equations :

$$\left( \frac{\partial \langle W_h \rangle}{\partial t}, \Phi_i \right) + \left( \nabla \cdot \underline{F}(\langle W_h \rangle), \Phi_i \right) = - \left( \underline{\tau}^{\text{RANS}}(\langle W_h \rangle), \Phi_i \right).$$

A natural hybridation writes :

$$\left( \frac{\partial W_h}{\partial t}, \Phi_i \right) + \left( \nabla \cdot \underline{F}(W_h), \Phi_i \right) = -\theta \left( \underline{\tau}^{\text{RANS}}(\langle W_h \rangle), \Phi_i \right) - (1-\theta) \left( \underline{\tau}^{\text{LES}}(W'_h), \Phi'_i \right)$$

where $W_h$ denotes now the hybrid variables and $\theta$ is the blending function which varies between 0 and 1 and is defined by :

$$\theta = 1 - f_d(1 - \tanh(\xi^2))$$

$$\text{with} \quad \xi = \frac{\Delta}{l_{\text{RANS}}} \quad \text{or} \quad \xi = \frac{\mu_{\text{SGS}}}{\mu_{\text{RANS}}}$$

the shielding function $f_d$ being defined as in Subsection 2.2.

The blending function $\theta$ allows for a progressive switch from RANS to LES (DVMS in our case) where the grid resolution becomes fine enough to resolve a significant part of the local turbulence scales or fluctuations, i.e. computational regions suitable for LES-like simulations. Additionally, this blending function, thanks to the shielding function $f_d$, prevents the activation of the LES
mode in the boundary layer.

By way of example, the isocontours of the blending function $\theta$ for the flow past a circular cylinder are shown in Fig. 3. We can see in particular that in the boundary layer the RANS model is activated, while in the wake the LES approach is recovered. However, the near-wake is RANS whereas it should be LES. This indicates that the blending function needs further improvement or that the near-wake mesh is too coarse for LES computation.

<figure><img src="image_6.png" /><figcaption>FIGURE 3 – Isocontours of the blending function for the circular cylinder benchmark (Reynolds 10⁶): for $\theta = 1$ (in red) the RANS model is activated, wherever $0 < \theta < 1$ additional resolved fluctuations are computed, and in the limit $\theta \to 0$ (in blue) the LES approach (DVMS in this work) is recovered.</figcaption></figure>

## 2.6 Hybrid DDES/DVMS

The key idea of the proposed hybrid DDES/DVMS model [16] is to use the DVMS approach instead of the LES component of DDES in locations where this component is expected to be activated, especially in wake regions where the DVMS approach allows more accurate prediction.

Assuming that the semi-discretization of the DDES equations writes :

$$ \left( \frac{\partial W_h}{\partial t}, \Phi_i \right) + \left( \nabla \cdot F(\underline{W}_h), \Phi_i \right) = - \left( \underline{\tau}^{\text{DDES}}(\underline{W}_h), \Phi_i \right), $$

the hybrid equations are then defined by :

$$ \left( \frac{\partial \underline{W}_h}{\partial t}, \Phi_i \right) + \left( \nabla \cdot \underline{F}(\underline{W}_h), \Phi_i \right) = -\theta \left( \underline{\tau}^{\text{DDES}}(\underline{W}_h), \Phi_i \right) - (1-\theta) \left( \underline{\tau}^{\text{LES}}(\underline{W}_h'), \Phi_i' \right) $$

where $\underline{W}_h$ denotes the hybrid variables and $\theta$ is the blending function defined as in the previous subsection.

# 3 Numerical discretization

## 3.1 Numerical scheme

The governing equations are discretized in space using a mixed finite-volume/finite-element method applied to unstructured tetrahedral grids. The adopted scheme is vertex-centered, i.e.
all degrees of freedom are located at the vertices. The diffusive terms are discretized using P1 Galerkin finite-elements on the tetrahedra, whereas finite-volumes are used for the convective terms. The numerical approximation of the convective fluxes at the interface of neighboring cells is based on the Roe scheme [28] with low-Mach preconditioning [12]. To obtain second-order accuracy in space, the Monotone Upwind Scheme for Conservation Laws reconstruction method (MUSCL) [37] is used, in which the Roe flux is expressed as a function of reconstructed values of W at each side of the interface between two cells. A particular attention has been paid to the dissipative properties of the resulting scheme, since this is a key point for its successful use in LES, and therefore in simulations performed with a hybrid turbulence model. The numerical dissipation provided by the scheme used in the present work is made of sixth-order space derivatives [4] and thus is concentrated on a narrow-band of the highest resolved frequencies. This is expected to limit undesirable damping of the large scales by numerical dissipation. Moreover, a parameter $\gamma$ directly controls the amount of introduced numerical viscosity and can be explicitly tuned in order to reduce it to the minimal amount needed to stabilize the simulation. Time advancing is carried out through an implicit linearized method, based on a second-order accurate backward difference scheme and on a first-order approximation of the Jacobian matrix [22]. The resulting numerical discretization is second-order accurate both in time and space. It should be noted that the spatial discretization used in this work leads to a superconvergent approximation, i.e. the accuracy can be well above second-order for some Cartesian meshes. One can also add that the objective is not a higher-order convergence but a strong reduction of dissipation and a certain reduction of the dispersion in the general case of a non-Cartesian (but not too irregular) mesh.
## 3.2 Mesh adaptation for rotating machines

Our numerical model has been extended to rotor-stator simulation with a Chimera technique.

The Chimera method aims at solving partial differential equations by decomposition into subdomains with overlap in order to avoid having to use a global mesh. This method allows the communication between the computational subdomains thanks to the overlapping of the subdomains. In our case we consider a decomposition in two domains as shown in Fig. 4, a fixed domain in red and a rotating domain in blue. For computations we start by locating the boundary nodes of domain 1 in domain 2, and reciprocally we locate the boundary nodes of domain 2 in domain 1, then the aerodynamic values of each boundary node are determined by interpolation and each domain performs its calculation with the new interpolated values.

<figure><img src="image_7.png" /><figcaption>FIGURE 4 – Definition of rotor and stator from an initial mesh : in red, the stator, in blue the rotor, and in gray, the overlap region.</figcaption></figure>

A mesh adaptation loop based on a metric depending of the Mach number has been combined with the Chimera algorithm. The Transient Fixed Point mesh adaptation algorithm [1, 2, 11] is applied to each domain from a metric field evaluated on the whole domain.
# 4 Applications

## 4.1 Flow past a cylinder

The predictions of the flow around a circular cylinder are presented. Two Reynolds numbers, $10^6$ and $2 \times 10^6$, based on the cylinder diameter, D, and on the freestream velocity, are considered. Only a few numerical investigations have been performed for Reynolds numbers higher than $10^6$. This interval is inside the supercritical regime which appears at Reynolds number higher than $2 \times 10^5$ and for which the separation becomes turbulent [34].

The computational domain is such that $-15 \le x, y \le 15$, and $-1 \le z \le 1$ where x, y and z denote the streamwise, transverse and spanwise directions respectively, the cylinder axis being located at $x = y = z = 0$. The mesh contains 4.8 millions nodes.

In order to control the computational costs, a wall law is applied in the close vicinity of the wall. For accuracy purpose, the Reichardt analytical law [13], which gives a smooth matching between linear, buffer and logarithmic regions, is chosen. As the $y^+$ normalized distance is generally subject to large variations in complex flows, the wall law is combined with low Reynolds modeling which locally damps the fully-turbulent model in regions in which the wall law does not cover the buffer zone.

<figure><img src="image_8.png" /><figcaption>FIGURE 5 – Computational domain</figcaption></figure>
Reynolds number 10⁶

The behavior of the hybrid models presented in Section 2 are first investigated in terms of flow bulk coefficients. From Table ??, it can be noted that the predictions of these coefficients are globally in good agreement with the experimental data and the numerical results in the literature. The lift time fluctuations are nevertheless better predicted by the RANS/DVMS model compared to the other two hybrid models. This good behavior is confirmed by the correct prediction of the distribution of the mean pressure coefficient, see Fig. 6. On the other hand, it can be observed from Table ?? and Fig. 6 that Smagorinsky and WALE SGS models give globally comparable results. It is also worth noting that the RANS model predicts rather well the bulk coefficients. Fig. 7 shows instantaneous isocontours of the vorticity magnitude for each of the turbulence models used for this benchmark. As one might expect, the unsteady RANS model is more dissipative and captures less flow detail than the hybrid models, with in particular a much more damped and regular wake.

<table><caption>TABLE 1 – Bulk coefficient of the flow around a circular cylinder at Reynolds number 1M, $\bar{C}_d$ holds for the mean drag coefficient, $C'_l$ is the root mean square of lift time fluctuation, $\bar{C}_{pb}$ is the pressure coefficient at cylinder basis, $L_r$ is the mean recirculation length, $\bar{\theta}$ is the mean separation angle.</caption><thead><tr><th></th><th>Mesh size</th><th>y<sub>w</sub><sup>+</sup></th><th>y<sub>m</sub><sup>+</sup></th><th>C&#773;<sub>d</sub></th><th>C'<sub>l</sub></th><th>-C&#773;<sub>pb</sub></th><th>L<sub>r</sub></th><th>&theta;&#773;</th><th>St</th></tr></thead><tbody><tr><td colspan="10"><strong>Present simulation</strong></td></tr><tr><td>RANS k &minus; &epsilon;</td><td>4.8M</td><td>20</td><td>100</td><td>0.23</td><td>0.06</td><td>0.23</td><td>0.69</td><td>132</td><td>0.46</td></tr><tr><td>DDES k &minus; &epsilon; Goldberg</td><td>4.8M</td><td>20</td><td>100</td><td>0.20</td><td>0.04</td><td>0.22</td><td>0.87</td><td>138</td><td>0.13</td></tr><tr><td colspan="10"><strong>DDES/DVMS</strong></td></tr><tr><td>k &minus; &epsilon; / cubic Smagorinsky</td><td>4.8M</td><td>20</td><td>100</td><td>0.20</td><td>0.02</td><td>0.22</td><td>0.82</td><td>135</td><td>0.42</td></tr><tr><td>k &minus; &epsilon; / cubic WALE</td><td>4.8M</td><td>1</td><td>100</td><td>0.20</td><td>0.02</td><td>0.26</td><td>0.80</td><td>132</td><td>0.58</td></tr><tr><td colspan="10"><strong>RANS / DVMS</strong></td></tr><tr><td>k &minus; &epsilon; / cubic Smagorinsky</td><td>4.8M</td><td>20</td><td>100</td><td>0.24</td><td>0.05</td><td>0.22</td><td>0.62</td><td>133</td><td>0.42</td></tr><tr><td>k &minus; &epsilon; / cubic Smagorinsky</td><td>4.8M</td><td>1</td><td>100</td><td>0.25</td><td>0.09</td><td>0.25</td><td>0.64</td><td>132</td><td>0.46</td></tr><tr><td>k &minus; &epsilon; / cubic WALE</td><td>4.8M</td><td>1</td><td>100</td><td>0.26</td><td>0.11</td><td>0.22</td><td>0.65</td><td>134</td><td>0.42</td></tr><tr><td colspan="10"><strong>Other simulations</strong></td></tr><tr><td>SRANS Catalano [5]</td><td>2.3M</td><td>-</td><td>-</td><td>0.39</td><td>-</td><td>0.33</td><td>-</td><td>-</td><td>-</td></tr><tr><td>SLES Catalano [5]</td><td>2.3M</td><td>-</td><td>-</td><td>0.31</td><td>-</td><td>0.32</td><td>-</td><td>-</td><td>-</td></tr><tr><td>URANS Catalano [5]</td><td>2.3M</td><td>-</td><td>-</td><td>0.40</td><td>-</td><td>0.41</td><td>1.37</td><td>-</td><td>0.31</td></tr><tr><td>LES Kim [18]</td><td>6.8M</td><td>-</td><td>-</td><td>0.27</td><td>0.12</td><td>0.28</td><td>-</td><td>108</td><td>-</td></tr><tr><td colspan="10"><strong>Measurements</strong></td></tr><tr><td>Exp. [30]</td><td></td><td></td><td></td><td>0.24</td><td>-</td><td>0.33</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Exp. [29]</td><td></td><td></td><td></td><td>0.22</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Exp. [35]</td><td></td><td></td><td></td><td>0.25</td><td>-</td><td>0.32</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Exp. [9]</td><td></td><td></td><td></td><td>-</td><td>-</td><td>-</td><td>-</td><td>130</td><td>-</td></tr><tr><td>Exp. [38]</td><td></td><td></td><td></td><td>0.2-0.4</td><td>0.1-0.15</td><td>0.2-0.34</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>
<figure><img src="image_9.png" /><figcaption>FIGURE 6 – Distribution over the cylinder surface of the mean pressure coefficient. Comparison between experimental data and numerical results at Reynolds $10^6$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 7 – Circular cylinder, Reynolds $10^6$: instantaneous isocontours of the vorticity magnitude. From top to bottom: RANS, DDES, DDES/DVMS and RANS/DVMS.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>FIGURE 8 – Lift and drag coefficient fluctuation. From top to bottom : RANS, DDES, DDES/DVMS and RANS/DVMS.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>FIGURE 9 – Instantaneous vorticity magnitude. From top to bottom : DDES, DDES/DVMS and RANS/DVMS.</figcaption></figure>
Reynolds number $2 \times 10^6$

The main outputs obtained by the RANS, DDES, DDES/DVMS and RANS/DVMS models are summed up in Table 2. Regarding the two hybrid models using the DVMS approach, only the results obtained with the Smagorinsky SGS model are shown, the WALE SGS model giving very similar results. It can be noticed that all models overpredict the separation angle, and that the fluctuations of the lift coefficient are overpredicted by the RANS and DDES/DVMS model. The RANS/DVMS approach gives on the whole the most satisfactory results. In Fig. 10, which shows the distribution of the mean pressure coefficient, the numerical results obtained with both hybrid models RANS/DVMS and DDES/DVMS are in very good agreement with the experimental results, in particular RANS/DVMS.

<table><caption>TABLE 2 – Bulk coefficient of the flow around a circular cylinder at Reynolds number $2 \times 10^6$. The subscript *S* holds for Smagorinsky SGS model.</caption><thead><tr><th></th><th>Mesh size</th><th>y<sub>w</sub><sup>+</sup></th><th>y<sub>m</sub><sup>+</sup></th><th>C&#773;<sub>d</sub></th><th>C'<sub>l</sub></th><th>-C&#773;<sub>pb</sub></th><th>L<sub>r</sub></th><th>&theta;&#773;</th><th>St</th></tr></thead><tbody><tr><td colspan="10"><strong>Present simulation</strong></td></tr><tr><td>RANS k &minus; &epsilon;</td><td>4.8M</td><td>20</td><td>100</td><td>0.26</td><td>0.066</td><td>0.30</td><td>-</td><td>128</td><td>-</td></tr><tr><td>DDES k &minus; &epsilon; Goldberg</td><td>4.8M</td><td>20</td><td>100</td><td>0.28</td><td>0.038</td><td>0.27</td><td>-</td><td>132</td><td>-</td></tr><tr><td colspan="10"><strong>DDES/DVMS</strong></td></tr><tr><td>k &minus; &epsilon; / cubic Smagorinsky</td><td>4.8M</td><td>40</td><td>100</td><td>0.26</td><td>0.026</td><td>0.35</td><td>0.83</td><td>130</td><td>0.33</td></tr><tr><td>k &minus; &epsilon; / cubic WALE</td><td>4.8M</td><td>40</td><td>100</td><td>0.24</td><td>0.020</td><td>0.30</td><td>1.15</td><td>128</td><td>0.19</td></tr><tr><td colspan="10"><strong>RANS/DVMS</strong></td></tr><tr><td>k &minus; &epsilon; / cubic Smagorinsky</td><td>4.8M</td><td>40</td><td>100</td><td>0.24</td><td>0.030</td><td>0.30</td><td>0.80</td><td>132</td><td>0.53</td></tr><tr><td>k &minus; &epsilon; / cubic WALE</td><td>4.8M</td><td>40</td><td>100</td><td>0.26</td><td>0.057</td><td>0.30</td><td>0.75</td><td>128</td><td>0.46</td></tr><tr><td colspan="10"><strong>Other simul.</strong></td></tr><tr><td>LES/TBLE [33]</td><td></td><td></td><td></td><td>0.24</td><td>0.029</td><td>0.36</td><td>-</td><td>105</td><td></td></tr><tr><td colspan="10"><strong>Measurements</strong></td></tr><tr><td>Exp. [30]</td><td></td><td></td><td></td><td>0.26</td><td>0.033</td><td>0.40</td><td>-</td><td>105</td><td></td></tr><tr><td>Exp. [29]</td><td></td><td></td><td></td><td>0.32</td><td>0.029</td><td>-</td><td>-</td><td>-</td><td></td></tr></tbody></table>
<figure><img src="image_13.png" /><figcaption>FIGURE 10 – Distribution over the cylinder surface of the mean pressure coefficient. Comparison between experimental data and numerical results at Reynolds 2 × 10⁶.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 11 – Lift and drag coefficient fluctuation. From top to bottom : DDES/DVMS and RANS/DVMS.</figcaption></figure>
## 4.2 Flow around a rotating cross

We present a preliminary computation obtained with the combination of the mesh adaptation, the Chimera approximation and the DDES model. They are used for simulating the mixing process obtained with the rotation of a cross in a cylindric box. The Reynolds number is 1.8 million, based on the thickness of the blades. The number of vertices of the initial mesh is 10.466. After two turns, and five mesh adaptations, the number of vertices of the adapted mesh is 1.05 million.

<figure><img src="image_15.png" /><figcaption>FIGURE 12 – DDES computation of a rotating mixing device with a mesh adaptive Chimera approach (horizontal cut). View of the mesh and of Q-factor (colored with velocity magnitude) after the vertical upper blade.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIGURE 13 – DDES computation of a rotating mixing device with a mesh adaptive Chimera approach (horizontal cut). Partial view of the mesh and velocity magnitude.</figcaption></figure>

In Fig. 12a, we present a view of this final adapted mesh. Mesh size normal to the cylinder is yet not small enough to allow the creation of strong non-2D features. This is verified from the examination of the Q-factor, Fig. 12b. Fig. 13 shows details of mesh and velocity magnitude. It seems from this first calculation that despite the high Reynolds number, the flow separates at blade angle is a rather stable mode.

## 5 Conclusion

Several hybrid strategies, based on the DDES, RANS and DVMS models, are first evaluated for the simulation of the supercritical flow around a circular cylinder. This benchmark is characterized by turbulent boundary-layer separation. The cylinder flow is first computed at $Reynolds\ 10^6$. The predictions of the main flow parameters are in overall good agreement with experimental data and the results of simulations in the literature, especially the RANS/DVMS model. A second series of computations of the flow past a circular cylinder is carried out at $Reynolds\ 2 \times 10^6$.
The predictions obtained by the hybrid models used to perform this benchmark are a little less accurate compared to the previous Reynolds number, more specifically for the mean separation angle which is overestimated. Again, the RANS/DVMS model is the hybrid model among those applied which gives the best results overall. It is also to be noted that the pressure distribution is well predicted by the hybrid models for both Reynolds numbers, especially by the RANS/DVMS model. Then, the DDES formulation is applied to the flow around a rotating cross inside a cylinder with a Reynolds number of $1.8 \times 10^6$. Preliminary results using mesh adaptation with meshes of 1 million vertices show that separations are more stable than in the flows around the cylinder. RANS/DVMS computations are planned. The new tool is also being applied to other rotating devices like propellers and will produce mesh adapted results for these geometries in a near future.

## Acknowledgements

This work was supported by the ANR NORMA project, grant ANR-19-CE40-0020-01 of the French National Research Agency. The authors gratefully acknowledge GENCI for granted access to HPC resources through CINES (grants 2020-A0092A05067 and 2021-A0102A06386) and IDRIS (grant 2021-A0112A05067).

## Références

[1] F. Alauzet, P.J. Frey, P.-L. George, and B. Mohammadi. 3D transient fixed point mesh adaptation for time-dependent problems : Application to CFD simulations. *J. Comp. Phys.*, 222 :592-623, 2007.

[2] F. Alauzet, A. Loseille, and G. Olivier. Multi-scale anisotropic mesh adaptation for time-dependent problems. RR-8929, INRIA, June 2016.

[3] P.A. Bakhvalov and T.K. Kozubskaya. On using artificial viscosity in edge-based schemes on unstructured grids. *Mathematical Models and Computer Simulations*, 13(4) :705-715, 2021.

[4] S. Camarri, M. V. Salvetti, B. Koobus, and A. Dervieux. A low diffusion muscl scheme for les on unstructured grids. *Comput Fluids*, 33 :1101-1129, 2004.

[5] P. Catalano, Meng Wang, G. Iaccarino, and P. Moin. Numerical simulation of the flow around a circular cylinder at high Reynolds numbers. *Int. J. Heat Fluid Flow*, 24 :463-469, 2003.

[6] A.P. Duben, T.K. Kozubskaya, P.V. Rodionov, and V.O. Tsvetkova. Ebr schemes with curvilinear reconstructions of variables in the near-wall region. *Computational Mathematics and mathematical physics*, 61(1) :1-16, 2021.

[7] V.A. Garanzha and L.N. Kudryavtseva. Moving deforming mesh generation based on quasi-isometric functional. *Lecture Notes in Computational Science and Engineering, WoS, Scopus-Q1*, 2020.

[8] M. Germano, U. Piomelli, P. Moin, and W. Cabot. A dynamic subgrid-scale eddy viscosity model. *Phys. Fluids A*, 3(7) :1760-1765, 1991.

[9] B. Goelling. Experimental investigations of separating boundary-layer flow from circular cylinder at Reynolds numbers from $10^5$ up to $10^7$; three-dimensional vortex flow of a circular cylinder. In G.E.A. Meier and K.R. Sreenivasan, editors, *Proceedings of IUTAM Symposium on One Hundred Years of Boundary Layer Research*, pages 455-462, The Netherlands, 2006. Springer.
[10] U. Goldberg, O. Peroomian, and S. Chakravarthy. A wall-distance-free $k-\epsilon$ model with enhanced near-wall treatment. *Journal of Fluids Engineering*, 120:457-462, 1998.

[11] D. Guégan, O. Allain, A. Dervieux, and F. Alauzet. An $L^\infty-L^p$ mesh adaptive method for computing unsteady bi-fluid flows. *Int. J. Numer. Meth. Eng.*, 84(11):1376-1406, 2010.

[12] H. Guillard and C. Viozat. On the behaviour of upwind schemes in the low Mach number limit. *Comput Fluids*, 28:63-86, 1999.

[13] J. Hinze. *Turbulence*. McGraw-Hill, 1959.

[14] Jens Holmen, Thomas J. R. Hughes, Assad A. Oberai, and Garth N. Wells. Sensitivity of the scale partition for variational multiscale large-eddy simulation of channel flow. *Physics of Fluids*, 16(3):824-827, 2004.

[15] T.J.R. Hughes, L. Mazzei, and K.E. Jansen. Large eddy simulation and the variational multiscale method. *Comput. Vis. Sci.*, 3:47-59, 2000.

[16] E. Itam, S. Wornom, B. Koobus, and A. Dervieux. Combining a ddes model with a dynamic variational multiscale formulation. In *12th International ERCOFTAC Symposium on Engineering Turbulence Modelling and Measurements (ETMM12)*, La Grande Motte, France, 2018.

[17] Marc Jaeger and Gouri Dhatt. An extended k-epsilon finite element model. *International Journal for Numerical Methods in Fluids*, 14(11):1325-1345, June 1992.

[18] S.-E. Kim and L.S. Mohan. Prediction of unsteady loading on a circular cylinder in high reynolds number flows using large eddy simulation. *Proceedings of OMAE 2005: 24th International Conference on Offshore Mechanics and Artic Engineering, June 12-16, Halkidiki, Greece, OMAE 2005-67044*, 2005.

[19] B. Koobus and C.. Farhat. A variational multiscale method for the large eddy simulation of compressible turbulent flows on unstructured meshes-application to vortex shedding. *Comput. Methods Appl. Mech. Eng.*, 193:1367-1383, 2004.

[20] Marie-Hélène Lallemand, Hervé Stève, and Alain Dervieux. Unstructured multigridding by volume agglomeration : Current status. *Computers & Fluids*, 21:397-433, 1992.

[21] D.K. Lilly. A proposed modification of the Germano subgrid-scale closure method. *Phys. Fluids*, A4:633, 1992.

[22] R. Martin and H. Guillard. A second-order defect correction scheme for unsteady problems. *Comput Fluids*, 25(1):9-27, 1996.

[23] F. Miralles, B. Sauvage, S. Wornom, B. Koobus, and A. Dervieux. Application of hybrid RANS/VMS modeling to rotating machines. In *Conference on Modelling Fluid Flow (CMFF'22), The 18th International Conference on Fluid Flow Technologies, Budapest, Hungary*, August 30-September 2, 2022.

[24] C. Moussaed, M. V. Salvetti, S. Wornom, B. Koobus, and A. Dervieux. Simulation of the flow past a circular cylinder in the supercritical regime by blending rans and variational-multiscale les models. *Journal of Fluids and Structures*, 47:114-123, 2014.

[25] C. Moussaed, S. Wornom, M.V. Salvetti, B. Koobus, and A. Dervieux. Impact of dynamic subgrid-scale modeling in variational multiscale large-eddy simulation of bluff body flows. *Acta Mechanica*, 225:3309-3323, 2014.

[26] F. Nicoud and F. Ducros. Subgrid-scale stress modelling based on the square of the velocity gradient tensor. *Flow, Turbulence and Combustion*, 62:183-200, 1999.

[27] H. Reichardt. Vollständige darstellung der turbulenten geschwindigkeitsverteilung in glatten leitungen. *Zeitschrift für Angewandte Mathematik und Mechanik*, 31:208-219, 1951.
[28] P.L. Roe. Approximate Riemann solvers parameters vectors and difference schemes. *J Comput Phys*, 43:357-372, 1981.

[29] G. Schewe. On the force fluctuations acting on a circular cylinder in crossflow from subcritical up to transcritical Reynolds numbers. *Journal of Fluid Mechanics*, 133:265-285, 1995.

[30] W.C.L. Shih, C. Wang, D. Coles, and A. Roshko. Experiments on flow past rough circular cylinders at large Reynolds numbers. *J. Wind Eng. Indust. Aerodyn.*, 49:351-368, 1993.

[31] J. Smagorinsky. General circulation experiments with the primitive equations. *Monthly Weather Review*, 91(3):99-164, 1963.

[32] P. R. Spalart, S. Deck, M. L. Shur, K. D. Squires, Strelets M. Kh., and A. K. Travin. A new version of detached-eddy simulation, resistant to ambiguous grid densities. *Theoritical and Computational Fluid Dynamic*, 20(3):181-195, 2006.

[33] A. Sreenivasan and B. Kannan. Enhanced wall turbulence model for flow over cylinder at high reynolds number. *AIP Advances*, 095012, 2019.

[34] B.M. Sumer and J. Fredsoe. *Hydrodynamics Around Cylindrical Structures*. World Scientific, 2006.

[35] E. Szechenyi. Supercritical Reynolds number simulation for two-dimensional flow over circular cylinders. *J. Fluid Mech.*, 70:529-542, 1975.

[36] Hubert Baya Toda, Karine Trun, and Franck Nicoud. Is the dynamic procedure appropriate for all sgs models. 2010.

[37] B. Van Leer. Towards the ultimate conservative scheme. iv : a new approach to numerical convection. *J Comput Phys*, 23:276-299, 1977.

[38] M. M. Zdravkovich. Different modes on vortex shedding : an overview. *Journal of Fluid and Structures*, 10(5):427-437, 1996.
informatics mathematics
Inria

RESEARCH CENTRE
SOPHIA ANTIPOLIS – MÉDITERRANÉE

2004 route des Lucioles - BP 93
06902 Sophia Antipolis Cedex

Publisher
Inria
Domaine de Voluceau - Rocquencourt
BP 105 - 78153 Le Chesnay Cedex
inria.fr

ISSN 0249-6399