Bruno Jeanneret
email: bruno.jeanneret@ifsttar.fr

Daniel Ndiaye

Sylvain Gillet

Rochdi Trigui

H&HIL: A novel tool to test control strategy with Human and Hardware In the Loop

With this work, the authors try to make HIL simulation more realistic with the introduction of the Human driver in the loop. To succeed in this objective we develop a set of tools to connect easily a wide variety of real or virtual devices together : of course the driver but also racing gamer joystick, real engine, virtual drivetrain, virtual driver environment ... The detailed approach represent a step forward before testing control strategy or new powertrain on a real vehicle. First results showed a good effectiveness and modularity of the tool.

I. INTRODUCTION

The transportation sector is responsible for a wide share of Energy consumption and pollutant emission everywhere in the world. The growth of environmental awareness is today among the most constringent drivers that the researchers and the manufacturers have to consider when designing environmental friendly solutions for drive trains, including electric, hybrid and fuel cell Vehicles. Nevertheless, most of the development and the evaluation of the new solutions still follow classical scheme that includes modeling and testing components and drivetrain under standard driving cycles or in the best case, using preregistered real world driving cycle [START_REF] Oh | Evaluation of motor characteristics for hybrid electric vehicles using the hil concept[END_REF], [START_REF] Trigui | Performance comparison of three storage systems for mild hevs using phil simulation[END_REF], [START_REF] Bouscayrol | Hardware in the loop simulation in "Control and Mechatronics[END_REF], [START_REF] Shidore | Phev 'all electric range' and fuel economy in charge sustaining mode for low soc operation of the jcs vl41m li-ion battery using battery hil[END_REF], [START_REF] Castaings | Comparison of energy management strategies of a battery/supercapacitors system for electric vehicle under real-time constraints[END_REF] and [START_REF] Verhille | Hardwarein-the-loop simulation of the traction system of an automatic subway[END_REF]. This methodology of design does not include the variability of the driving conditions generated by the two major factors that are the driver behavior and the infrastructure influence (elevation, turns, speed limitation, traffic jam, traffic lights, . . . ).

Moreover, the progress registered in the Intelligent Transportation Systems (ITS) makes it possible today to optimize the actual use of the vehicles (classical, EVs, HEVs) by capturing instantaneous information about the infrastructure (road profile, road signs), the traffic and also about the weather conditions. These informations are used for HEV and PHEVs to optimize their energy management and for mono-sources vehicles (classical, EV) to give inputs to the ADAS systems in order to perform a lower consumption and emission (ecodriving concepts for example). The simulation and testing using standard or preregistered driving cycles is therefore too limited to take into account all these aspects. The need is then to develop new simulation and testing schemes able to consider a more realistic modeling of the vehicle environment and to include the driver or its model in the simulation/testing loop.

The methodology presented in this paper is based on a progressive and modular approach for simulating and testing in a HIL configuration different types of drive trains while including the vehicle environment models and the driver. The modularity is developed considering two axes:

• virtual to real : The ability of the developed system ranges from all simulation configuration (SIL) to driver + hardware in the loop simulation (H&HIL)

• multiplatform: The communication protocol between different systems or models allows easy exchange of systems and simulation platforms (energetic and dynamic vehicle model, different driving simulators configuration, Power plant in a HIL configuration, . . . ). The chosen protocol, namely the CAN (Control Area Network) can also easily be used to address properly further steps like vehicle prototype design

In the following sections we will describe the developed concepts and tools. The last chapter presents two applications made with the tool. A conclusion lists applications that can be developped with the facility.

II. TOOL DESCRIPTION

A. Architecture of the tool

The main program named MODYVES aims at connecting a generic input and a driver to a generic output with a vehicle model. The main architecture of the framework is presented in figure 1. Inputs can be choosen among:

• a Gamepad like Logitech G27 Joystick or equivalent. Modyves uses Simple DirectMedia Layer1 to provide low level access to the device

B. Communication layer

The CAN protocol is widely used in automotive industry. For this reason, we selected this kind of network for the communication of the different pieces of hardware. Easy to implement, another advantage is that the communication can be secured by testing the Rx time (Received frame) in a real time clock.

Two different USB to CAN converter have been integrated in the tool, PEAK2 module and Systec3 module.

Both provide device driver (Dynamic Link Library files) and header files to connect to the module and parametrise it, decode received frame and send transmitted ones.

C. Hardware

Two main hardwares have been integrated in our tool, a prototyping hardware commercialised by dSPACE, namely Micro-autobox (MABX-II), and a low cost 32 bit microcontroler from Texas Instrument, the C2000 Peripheral Experimenter kit equipped with a F28335 control card. Both are supported by The Mathworks with matlab and simulink coder dSPACE hardware, and matlab and embedded coder for C2000.

The table I makes a short comparison between hardwares. Of course these two devices have not the same performances and are not suitable for the same kind of application.

One can notice the power of dSPACE board (procesor speed and size of RAM and Flash), but the major advantage of the tool is its combinations with a RTI interface and controlDesk software that enable to rapidly develop, debug and validate a real time project.

On the other hand, TI C2000 control card is well suited for a variety of automotive applications. If offers a good processor speed with a moderate but sufficient Flash memory to developp and embed a real time application at a moderate price. Debugging is far more complex than with dSPACE product.

D. Software

Modyves is written in python and is cadenced by a windows timer, so depending upon computer characteristics and computer load, it can deviate from its theorical execution period. Nevertheless, when connected to real plant, the only critical part running in Modyves is the driver behavior and the communication layer.

The IFSTTAR driving simulator is a static driving simulator based on a real Peugeot 308 and a software part named SIM2. The cockpit and all commands are unchanged to offer the most realistic driving environments to the driver. An embedded electronic card in the vehicle reads all sensors values placed on the pedals, the gearbox and the steering wheel. It also controls the force feedback on the steering wheel. The electronic card controls the vehicle dashboard by sending CAN informations through the OBD connector. SIM2 is the IFSTTAR simulator software, and contains various types of models: vehicle, road, sound, traffic and visual models. The road scene is displayed on 5 screens offering up to 180°h

orizontal Field Of View. The IFSTTAR driving simulator is used in the field of human factor research, ergonomics studies, energy-efficient driving, advanced training and studies of "Human and Hardware in Loop"( H&HIL).

VEHLIB is namely a Simulink library. A framework has been developped for years around this library to integrate all the necessary components model to develop and simulate conventional, hybrid or all electric vehicles [START_REF] Vinot | Model simulation, validation and case study of the 2004 THS of Toyota Prius[END_REF] [START_REF] Jeanneret | New Hybrid Concept Simulation Tools, evaluation on the Toyota Prius Car[END_REF]. VEHLIB is developped with a combined backward or forward approach. The forward models are able to run on real time hardware with their respectives connexion blocks to real environment [START_REF] Jeanneret | Mise en oeuvre d une commande temps reel de transmission hybride sur banc moteur[END_REF].

III. EXAMPLES OF APPLICATION

The figure 3 presents different step of integration from the pure simulation on a personal computer up to the deployement on a final facilty. This procedure has been succesfully conducted for IFSTTAR's driving simulator. For simplicity reason, the vehicle model has been compiled in a C2000 hardware instead of linking the dynamic library to the simulator software. As presented in figure 1, the tool allows to make a large variety of tests and integration depending on the effort and the final objective. Two applications are presented hereafter, the first one is a model in the loop application with a G27 joystick, the second is a Power Hardware in the Loop with a driver simulator, a virtual vehicle connected to a real engine. Both include a driver in the loop.

A. G27 Joystick and Model in the loop application

A first example to introduce the human in the loop consists of a real time application running in a soft real time mode (i.e. a windows timer). It is easy to implement and need only a G27 Joystick and a computer to run (See example 1 in figure II-D).

As mentioned earlier, this application uses a windows timer to set the switching time of the application. In order to verify the real execution period of Modyves, the jitter time has been monitored on a personal laptop (intel core i7 3610QM @2.3 GHz running Windows 7) without perturbation (no other program were running on the computer). The theoretical frequency is 1 kHz. The results are presented in figure 4. The mean step time is 0.001002 sec., the max value is 0.0027 sec. (only one occurence during this test which lasted around 30 seconds) and the minimum value is 0.0009999 second. These values do not affect neither the driver perception, nor the model behavior, because they are far from the time constant of these "systems". Consequently the deviation from the theoretical period is small enough for our application.

B. Driver simulator with HIL application on engine test bench

This is the most complete situation described in figure II-D. In this case, the vehicle model (excepted the engine) is running in a hard real time platform, a Micro-autobox hardware in this particular case. The latter communicates with the engine test bench and exchange some informations with it, namely:

• send accelerator pedal position to the engine and rotation speed to the electric generator

• receive actual torque

The bench is presented in figure 7. The bench is running in the so-called throttle/speed mode. At each time step, the actual torque is measured on the bench, transmitted to the model and introduced in the simulated clutch. It passes throwgh the different components up to the vehicle wheels. The longitudinal motion equation is resolved allowing to calculate the speed of the different shafts up to the engine speed, which is in turn transmitted to the bench as a target generator speed. In the same time, the accelerator position is also send to the engine ECU.

The driving simulator is presented in figure 8. One can notice that fuel cut-off is effective in engine simulation when the vehicle decelerate (figure 6) but are not present in the bench (figure 10) where cut-off is not validated in the actual engine ECU.

IV. CONCLUSION

In this paper, a generic framework has been presented to develop and test mecatronic applications in a progressive way. The facility is used in the laboratory to test ADAS (Advanced Driver Assistance Systems) and allows to take into account driver behavior as well as realistic vehicle environment. In this context, not only the fuel consumption is measured in a real engine, but one can measure engine emissions thanks to the gas analyzers available in the laboratory including CO, HC, NOx and particles measurements.

Number of applications could be performed with this facility:

• the coupled setting could be used to approach a real use behavior of the vehicle and the ICE instead of performing only standard driving cycles. In fact, the new emission regulation consider pollutant measurement on a real track in a procedure called Real Drive Emission (RDE) using Portable devices (PEMS). In order to anticipate this phase, the coupled test bench could help for the engine design and tuning in order to reduce real emissions

• the facility could embed ADAS system. They can be easily implemented and tested on urban, road or highway context with a human driver in the loop. For example, the assistance systems for eco-driving could be assessed in terms of near-real fuel consumption and pollutant emission

• the simulator can be a platform to evaluate different degrees of driving delegation towards autonomous vehicle. For example, it is easy to simulate manual or automatic gearbox, to test different kind of speed regulator, speed limitation...

• hybrid vehicles can be emulated in the real time model or they can be introduced in the bench to simulate hybrid vehicle. A clutch can also be physically present to test all electric mode for a parallel hybrid vehicle. With this kind of application, several energy management laws can be quickly implemented and tested in a almost real environment A further step for the software can consist of generalizing it to allow the core module to connect to any kind of devices as described in figure 11.

Fig. 1 .

 1 Fig. 1. The Modyves framework

Fig. 2 .

 2 Fig. 2. VEHLIB in one H&HIL configuration

Fig. 3 .

 3 Fig. 3. The different step of integration

Fig. 4 .

 4 Fig. 4. Jitter for a 1 kHz application

Fig. 5 .Fig. 6 .

 56 Fig. 5. Vehicle speed and driver actuators

Fig. 7 .

 7 Fig. 7. Engine on the test bench

Fig. 8 .

 8 Fig. 8. Driving simulator Figures 9 and 10 illustrate the behaviour of the driver and the response of the engine on the test bench.

Fig. 9 .Fig. 10 .

 910 Fig. 9. Vehicle speed and driver actuators

Fig. 11 .

 11 Fig. 11. The engine for software integration

TABLE I .

 I HARDWARE CHARACTERISTICS

https://www.libsdl.org/

http://www.peak-system.com/

http://www.systec-electronic.com/


Bernd Kister

Stéphane Lambert1

Bernard Loup

IMPACT TESTS ON SMALL SCALE EMBANKMENTS WITH ROCKERY -LESSONS LEARNED

Keywords: rockfall protection embankment (RPE), impact, rockery, block shape, ratio rotational to translational energy, freeboard, energy dissipation

In the project AERES (Analysis of Existing Rockfall Embankments of Switzerland) smallscale quasi-2D-experiments had been done with embankments with stones placed parallel to the slope and stones placed horizontally. The experiments showed that rotating cylinders acting as impactors may surmount an embankment with batter 2:1, even if the freeboard is chosen to 1.5 times the block diameter. So a slope with an inclination of about 60° and equipped with rockery in general does not guarantee that a freeboard of one block diameter will be sufficient as described by the Austrian standard. During the test series also a block with an octagonal cross section had been used. This block, with no or only very low rotation, on the other hand was not able to surmount an embankment with rockery and a freeboard of about 0.8 times the block diameter. The evaluation of test data showed additional that the main part of energy dissipation occurs during the first 6 ms of the impact process. At least 75% to 85% of the block's total kinetic energy will be transformed into compression work, wave energy and heat when the block hits the embankment.

INTRODUCTION

A common measure used in Switzerland to stabilize steep slopes of rockfall protection embankments is the use of rockery. According to the interviews done with employees of canton departments as well as with design engineers the rockery used in the past at the uphill slope of an embankment was constructed with a batter between 60° to 80°. But in general only less attention had been paid on the behavior of such natural stone walls during the impact of a block. The main reasons to use rockery was to limit the area, which is necessary for embankment construction, and/or to stop rolling blocks.

To check the experimental results of [START_REF] Hofmann | Bemessungsvorschlag für Steinschlagschutzdämme[END_REF] and to study the impact process of blocks impacting rockfall protection embankments (RPE) with a rockery cover at the uphill slope, small scale quasi-2D-experiments have been executed and analyzed at the Lucerne University of Applied Sciences and Arts (HSLU) during the project AERES [START_REF] Kister | Analysis of Existing Rockfall Embankments of Switzerland (AERES)[END_REF].

The load case "impact of a block onto an embankment" may be divided into two scenarios:

 The embankment is punched-through by the block and  the embankment is surmounted by the block. The first one is a question of stability of the construction itself while the second one concerns the fitness for purpose of a RPE. The most significant results of the tests done in the project AERES concerning the fitness for purpose of a RPE are shown below.

TEST CONDITIONS

Two types of rockery with different orientation of stones at the "uphill" slope had been studied. To create a relatively "smooth" rockery surface the stones had been placed parallel to the slope (Fig. 1a). For a "rough" rockery surface the stones had been placed horizontally, which resulted in a graded slope (Fig. 1b). Additional impact tests had been done on an embankment with a bi-linear slope with rockery at the lower part and soil at the upper part (Fig. 1c). The batter of the "rockery" was chosen to be 2:1, 5:2 and 5:1 in the tests. Three types of impactors had been used: A concrete cylinder G with a ratio of rotational to translational energy > 0.3, a hollow cylinder GS with a triaxial acceleration sensor inside with a ratio of rotational to translational energy between 0.2 and 0.1 and a block OKT with an octagonal cross section with no or only very low rotation. The ratio of rotational to translational energy of block GS corresponds very well with the results of in-situ tests done with natural blocks [START_REF] Usiro | An experimental study related to rockfall movement mechanism[END_REF]. The impact translational velocity in most tests was between 6 m/s and 7 m/s. Transformed to a prototype embankment with height of approx. 7 m this will result in real world block velocities of about 18 m/s to 21 m/s [START_REF] Kister | Analysis of Existing Rockfall Embankments of Switzerland (AERES)[END_REF].

FREEBOARD

The statements of Hofmann & Mölk [START_REF] Hofmann | Bemessungsvorschlag für Steinschlagschutzdämme[END_REF] concerning the freeboard have been transferred to the Austrian technical guideline ONR 24810:2013 [2] and it is said that the freeboard for an embankment with riprap (resp. rockery) and a slope angle of 50° or more should be at least one block diameter. To determine the maximum climbing height of a block during an impact and to get information about the influence of the roughness of the rockery surface, two impact

a) b) c)
tests on embankment models with a batter 2:1, but with different "rockery roughness" had been done. For these tests the freeboard was chosen to be approximately 1.9 times the block diameter. This value is a little bit less than the value of 2 times the block diameter specified in [2] for the freeboard of pure soil embankments, but significantly larger than the minimum value specified for embankments with rockery. The impact point was at a level, where the embankment thickness is larger than three times the block diameter (Fig. 2) and therefore according to [START_REF] Kister | Development of basics for dimensioning rock fall protection embankments in experiment and theory[END_REF] there was no risk that the embankment will be punched through. For "rough" rockery surface the climbing height of block GS was 1.8 times the block diameter for the first impact and 1.55 times the block diameter for the second impact (Fig. 2). The "rough" surface of the rockery led to a larger climbing height for the block than the "smooth" surface, although the block velocities had been very similar. The first impact of the block resulted in a larger climbing height than the second impact for both surface roughness types.

The tests showed that a slope with an inclination of about 60° and equipped with rockery in general does not guarantee that a freeboard of approx. one block diameter will be sufficient as described by the Austrian standard.

BLOCK SHAPE AND BLOCK ROTATION

During the test series the block OKT with an octagonal cross section and with no or only very low rotation was not able to surmount an embankment with rockery, if a crest to block diameter ratio of approx. 1.1 was chosen, even so the freeboard was only about 0.8 times the block diameter. Fig. 3 shows the trajectories of the concrete cylinder G with rotation and the block OKT impacting an embankment with stones placed horizontally at the "uphill" slope, batter of rockery 5:2. The difference in the impact translational velocities of both blocks was about 7%, which is within the measurement error (G: 5.9 m/s, OKT: 6.3 m/s). So block shape and block rotation play a significant role in the impact process.

ENERGY DISSIPATION

The evaluation of test data received in the project AERES showed that the main part of energy dissipation occurs during the first 6 ms of the impact process. During this period the block translational velocity is reduced to less than the half value for all three types of impactors used in the project. Differences in loss of block velocity and block energy of the three used impactors mainly occur after the large drop of velocity and energy.

CONCLUSION

The following parameters had been found, which are dominating the impact process and led the embankment be either surmounted by a block or punched-through:  The total block energy,  the ratio of rotational to translational block energy,  the impact angle, which is a function of block trajectory and slope inclination,  the shape of the block,  the embankment's thickness at the impact point. These parameters are mainly responsible for the fitness for purpose of a rockfall protection embankment. The experiments have shown that there are some interactions between these parameters, which could not be solved in detail with the existing experimental set-up. Further research has to be done to determine the freeboard in case of blocks with a natural shape and with a ratio of rotational to translational energy between the limits 0.1 and 0.2.

Fig. 1

 1 Fig. 1 Orientation of stones at the "uphill" slope: a) upright, b) horizontal, c) upright, upper slope without stones and reduced slope angle

Fig. 2 :

 2 Fig. 2: Max. climbing height CH max of impactor GS, stones placed horizontally, freeboard FB = 1.9*2r: a) first impact: CH max approx. 1.8*2r, b) second impact: CH max approx. 1.55*2r, 2r = block diameter

Fig. 3 :

 3 Fig. 3: Influence of block shape and block rotation: Trajectories of cylinder G (a) and block OKT (b), embankment with stones placed horizontally at the "uphill" slope, batter of rockery 5:2

Lucerne University of Applied Sciences and Arts, Technikumstrasse

21, CH -6048 Horw, Switzerland, since 2017: kister -geotechnical engineering & research, Neckarsteinacher Str. 4 B, D -Neckargemünd, +49 6223 71363, kister-ger@t-online.de 2 irstea, 2 rue de la Papeterie -BP 76, 38402 Saint-Martin-d'Hères cedex, France, +33 (0)4 76 76 27 94, stephane.lambert@irstea.fr 3 Federal Office for the Environment (FOEN), 3003 Bern, Switzerland, Tel. +41 58 465 50 98, bernard.loup@bafu.admin.ch


Pranav Chandramouli
email: pranav.chandramouli@inria.fr

Dominique Heitz

Sylvain Laizet

Etienne Mémin

Coarse large-eddy simulations in a transitional wake flow with flow models under location uncertainty

The focus of this paper is to perform coarse-grid large eddy simulation (LES) using recently developed sub-grid scale (SGS) models of cylinder wake flow at Reynolds number (Re) of 3900. As we approach coarser resolutions, a drop in accuracy is noted for all LES models but more importantly, the numerical stability of classical models is called into question. The objective is to identify a statistically accurate, stable sub-grid scale (SGS) model for this transitional flow at a coarse resolution. The proposed new models under location uncertainty (MULU) are applied in a deterministic coarse LES context and the statistical results are compared with variants of the Smagorinsky model and various reference data-sets (both experimental and Direct Numerical Simulation (DNS)). MULU are shown to better estimate statistics for coarse resolution (at 0.46% the cost of a DNS) while being numerically stable. The performance of the MULU is studied through statistical comparisons, energy spectra, and sub-grid scale (SGS) contributions. The physics behind the MULU are characterised and explored using divergence and curl functions. The additional terms present (velocity bias) in the MULU are shown to improve model performance. The spanwise periodicity observed at low Reynolds is achieved at this moderate Reynolds number through the curl function, in coherence with the birth of streamwise vortices.

Large Eddy Simulation, Cylinder Wake Flow

Introduction

Cylinder wake flow has been studied extensively starting with the experimental works of Townsend (1; 2) to the numerical works of Kravchenko and others (3; 4; 5). The flow exhibits a strong dependance on Reynolds number Re = U D/ν, where U is the inflow velocity, D is the cylinder diameter and ν is the kinematic viscosity of the fluid. Beyond a critical Reynolds number Re ∼ 40 the wake becomes unstable leading to the well known von Kármán vortex street. The eddy formation remains laminar until gradually being replaced by turbulent vortex shedding at higher Re. The shear layers remain laminar until Re ∼ 400 beyond which the transition to turbulence takes place up to Re = 10 5 -this regime, referred to as the sub-critical regime, is the focus of this paper.

The transitional nature of the wake flow in the sub-critical regime, especially in the shear layers is a challenging problem for turbulence modelling and hence has attracted a lot of attention. The fragile stability of the shear layers leads to more or less delayed roll-up into von Kármán vortices and shorter or longer vortex formation regions. As a consequence significant discrepancies have been observed in near wake quantities both for numerical simulations [START_REF] Ma | Dynamics and lowdimensionality of a turbulent near wake[END_REF] and experiments [START_REF] Parnaudeau | Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900[END_REF].

Within the sub-critical regime, 3900 has been established as a benchmark Re. The study of [START_REF] Parnaudeau | Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900[END_REF] provides accurate experimental data-set showing good agreement with previous numerical studies contrary to early experimental datasets [START_REF] Beaudan | Numerical experiments on the flow past a circular cylinder at sub-critical Reynolds number[END_REF]. The early experiments of Lourenco and Shih (9) obtained a V-shaped mean streamwise velocity profile in the near wake contrary to the U-shaped profile obtained by [START_REF] Beaudan | Numerical experiments on the flow past a circular cylinder at sub-critical Reynolds number[END_REF]. The discrepancy was attributed to inaccuracies in the experiment -a fact confirmed by the studies of [START_REF] Mittal | Suitability of Upwind-Biased Finite Difference Schemes for Large-Eddy Simulation of Turbulent Flows[END_REF] and (4). Parnaudeau et al.'s [START_REF] Parnaudeau | Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900[END_REF] experimental database, which obtains the U-shaped mean profile in the near wake, is thus becoming useful for numerical validation studies. With increasing computation power, the LES data sets at Re = 3900 have been further augmented with DNS studies performed by [START_REF] Ma | Dynamics and lowdimensionality of a turbulent near wake[END_REF].

The transitional nature of the flow combined with the availability of validated experimental and numerical data-sets at Re = 3900 makes this an ideal flow for model development and comparison. The LES model parametrisation controls the turbulent dissipation. A good SGS model should ensure suitable dissipation mechanism. Standard Smagorinsky model [START_REF] Smagorinsky | General circulation experiments with the primitive equations[END_REF] based on an equilibrium between turbulence production and dissipation, has a tendency to overestimate dissipation in general [START_REF] Meyers | On the model coefficients for the standard and the variational multi-scale Smagorinsky model[END_REF]. In transitional flows, where the dissipation is weak, such a SGS model leads to laminar regimes, for example, in the shear layers for cylinder wake. Different modifications of the model have been proposed to correct this behaviour. As addressed by [START_REF] Meyers | On the model coefficients for the standard and the variational multi-scale Smagorinsky model[END_REF], who introduced relevant improvements, the model coefficients exhibit a strong dependency both on the ratio between the integral length scale and the LES filter width, and on the ratio between the LES filter width and the Kolmogorov scale. In this context of SGS models, coarse LES remains a challenging issue.

The motivation for coarse LES is dominated by the general interest towards reduced computational cost which could pave the way for performing higher Re simulations, sensitivity analyses, and Data Assimilation (DA) studies. DA has gathered a lot of focus recently with the works of ( 13), [START_REF] Gronskis | Inflow and initial conditions for direct numerical simulation based on adjoint data assimilation[END_REF], and [START_REF] Yang | Enhanced ensemble-based 4dvar scheme for data assimilation[END_REF] but still remains limited by computational requirement.

With the focus on coarse resolution, this study analyses the performance of LES models for transitional wake flow at Re 3900. The models under location uncertainty (16; 17) are analysed in depth for their performance at a coarse resolution and compared with classical models. The models are so called as the equations are derived assuming that the location of a fluid parcel is known only up to a random noise i.e. location uncertainty. Within this reformulation of the Navier-Stokes equations, the contributions of the subgrid scale random component is split into an inhomogeneous turbulent diffusion and a velocity bias which corrects the advection due to resolved velocity field. Such a scheme has been shown to perform well on the Taylor Green vortex flow [START_REF] Harouna | Stochastic representation of the Reynolds transport theorem: revisiting large-scale modeling[END_REF] at Reynolds number of 1600, 3000, and 5000. The new scheme was shown to outperform the established dynamic Smagorinky model especially at higher Re. However, this flow is associated to an almost isotropic turbulence and no comparison with data is possible (as it is a pure numerical flow). Here we wish to assess the model skills with respect to more complex situations (with laminar, transient and turbulent areas) and coarse resolution grids. We provide also a physical analysis of the solutions computed and compare them with classical LES schemes and experimental data. Although the models are applied to a specific Reynolds number, the nature of the flow generalises the applicability of the results to a wide range of Reynolds number from 10 3 -10 5 , i.e. up to the pivotal point where the transition into turbulence of the boundary layer starts at the wall of the cylinder. The goal is to show the ability of such new LES approaches for simulation at coarse resolution of a wake flow in the subcritical regime. Note that recently for the same flow configuration, [START_REF] Resseguier | Stochastic modelling and diffusion modes for proper orthogonal decomposition models and small-scale flow analysis[END_REF] have derived the MULU in a reduced-order form using Proper Orthogonal Decomposition (POD), successfully providing physical interpretations of the local corrective advection and diffusion terms. The authors showed that the near wake regions like the pivotal zone of the shear layers rolling into vortices are key players into the modelling of the action of small-scale unresolved flow on the resolved flow .

In the following, we will show that the MULU are able to capture, in the context of coarse simulation, the essential physical mechanisms of the transitional very near wake flow. This is due to the split of the SGS contribution into directional dissipation and velocity bias. The next section elaborates on the various SGS models analysed in this study followed by a section on the flow configuration and numerical methods used. A comparison of the elaborated models and the associated physics is provided in the results section. Finally, a section of concluding remarks follows.

Models under location uncertainty

General classical models such as Smagorinsky or Wall-Adaptive Local Eddy (WALE) viscosity model proceed through a deterministic approach towards modelling the SGS dissipation tensor. However, [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF] suggests a stochastic approach towards modelling the SGS contributions in the Navier-Stokes (NS) equation. Building a stochastic NS formulation can be achieved via various methods. The simplest way consists in considering an additional additive random forcing [START_REF] Bensoussan | Equations stochastiques du type Navier-Stokes[END_REF]. Other modelling considered the introduction of fluctuations in the subgrid models (21; 22). Also, in the wake of Kraichnan's work [START_REF] Kraichnan | The structure of isotropic turbulence at very high Reynolds numbers[END_REF] another choice consisted in closing the large-scale flow in the Fourier space from a Langevin equation (24; 25; 26). Lagrangian models based on Langevin equation in the physical space have been also successfully proposed for turbulent dispersion [START_REF] Sawford | Generalized random forcing in random-walk models of turbu[END_REF] or for probability density function (PDF) modelling of turbulent flows (28; 29; 30). These attractive models for particle based representation of turbulent flows are nevertheless not well suited to a large-scale Eulerian modelling.

In this work we rely on a different stochastic framework of the NS equation recently derived from the splitting of the Lagrangian velocity into a smooth component and a highly oscillating random velocity component (i.e. the uncertainty in the parcel location expressed as velocity) [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF]:

dX t dt = u(X t , t) + σ(X t , t) Ḃ
The first term, on the right-hand side represents the large-scale smooth velocity component, while the second term is the small-scale component. This latter term is a random field defined from a Brownian term function Ḃ = dB/dt and a diffusion tensor σ. The small-scale component is rapidly decorrelating at the resolved time scale with spatial correlations (which might be inhomogeneous and non stationary) fixed through the diffusion tensor. It is associated with a covariance tensor:

Q ij (x, y, t, t ) = E((σ(x, t)dB t )(σ(x, t)dB t ) T ) = c ij (x, y, t)δ(t -t )dt. ( 1
)
In the following the diagonal of the covariance tensor, termed here as the variance tensor, plays a central role; it is denoted as a(x) = c(x, t). This tensor is a (3 × 3) symmetric positive definite matrix with dimension of a viscosity in m 2 s -1 . With such a decomposition, the rate of change of a scalar within a material volume, is given through a stochastic representation of the Reynolds Transport Theorem (RTT) (16; 17). For an incompressible small-scale random component (∇•σ = 0) the RTT has the following expression: d

V(t) q = V(t) d t q + ∇•(q u)- 1 2 d i,j=1 ∂ x i (a ij ∂ x j q) dt+∇q •σdB t dx. ( 2
)
where the effective advection u is defined as:

u = u - 1 2 ∇ • a. (3)
The first term on the right-hand side represents the variation of quantity q with respect to time: d t q = q(x, t+dt)-q(x, t). It is similar to the temporal derivative. It is important here to quote that q is a non differentiable random function that depends among other things on the particles driven by the Brownian component and flowing through a given location. The second term on the right-hand side stands for the scalar transport by the largescale velocity. However, it can be noticed that this scalar advection is not purely a function of the large-scale velocity. Indeed, the large-scale velocity is here affected by the inhomogeneity of the small-scale component through a modified large-scale advection (henceforth termed as velocity bias u ), where the effect of the fluctuating component is taken into account via the small-scale velocity auto-correlations a = (σσ T ). A similar modification of the large-scale velocity was also suggested in random walks Langevin models by [START_REF] Macinnes | Stochastic particle dispersion modeling and the tracer-particle limit[END_REF] who studied various stochastic models for particle dispersion -they concluded that an artificially introduced bias velocity to counter particle drift was necessary to optimise the models for a given flow. In the framework of modelling under location uncertainty, this term appears automatically. The third term in the stochastic RTT corresponds to a diffusion contribution due to the small-scale components. This can be compared with the SGS dissipation term in LES Modelling. This dissipation term corresponds to a generalization of the classical SGS dissipation term, which ensues in the usual context from the Reynolds decomposition and the Boussinesq's eddy viscosity assumption. Here it figures the mixing effect exerted by the smallscale component on the large-scale component. Despite originating from a very different construction, in the following, for ease of reference, we keep designating this term as the SGS contribution. The final term in the equation is the direct scalar advection by the small-scale noise.

It should be noted that the RTT corresponds to the differential of the composition of two stochastic processes. The Ito formulae, which is restricted to deterministic functions of a stochastic process, does not apply here. An extended formulae know as Ito-Wentzell (or generalized Ito) formulae must be used instead [START_REF] Kunita | Stochastic Flows and Stochastic Differential Equations[END_REF].

Using the Stochastic RTT, the large-scale flow conservation equations can be derived (for the full derivation please refer to (16; 17)). The final conservation equations are presented below:

Mass conservation:

d t ρ + ∇ • (ρ u)dt + ∇ρ •σ dB t = 1 2 ∇ • (a∇ρ)dt, (4)
which simplifies to the following constraints for an incompressible fluid:

∇ • σ = 0, ∇ • u = 0, (5)
The first constraint maintains a divergence free small-scale velocity field, while the second imposes the same for the large smooth effective component.

We observe that the large-scale component, u, is allowed to be diverging, with a divergence given by ∇ • ∇ • a. As we shall see, this value is in practice quite low. This weak incompressibility constraint results in a modified pressure computation, which is numerically not difficult to handle. Imposing instead a stronger incompressibility constraint on u introduces an additional cumbersome constraint on the variance tensor (∇ • ∇ • a = 0). In this work we will rely on the weak form of the incompressibility constraint. The largescale momentum equation boils down to a large-scale deterministic equation after separation between the bounded variation terms (i.e. terms in "dt") and the Brownian terms, which is rigorously authorized -due to uniqueness of this decomposition. Momentum conservation:

∂ t u + u∇ T (u - 1 2 ∇ • a) - 1 2 ij ∂ x i (a ij ∂ x j u) ρ = ρg -∇p + µ∆u. ( 6
)
Similar to the deterministic version of the NS equation, we have the flow material derivative, the forces, and viscous dissipation. The difference lies in the modification of the advection which includes the velocity bias and the presence of the dissipation term which can be compared with the SGS term present in the filtered NS equation. Both the additional terms present in the stochastic version are computed via the auto-correlation tensor a. Thus to perform a LES, one needs to model, either directly or through the smallscale noise, the auto-correlation tensor. Two methodologies can be envisaged towards this: the first would be to model the stochastic small-scale noise (σ(X t , t) Ḃ) and thus evaluate the auto-correlation tensor. We term such an approach as purely 'stochastic LES'. The second method would be to model the auto-correlation tensor directly as it encompasses the total contribution of the small scales. This method can be viewed as a form of 'deterministic LES' using stochastically derived conservation equations and this is the approach followed in this paper. The crux of the 'deterministic LES' approach thus revolves around the characterisation of the auto-correlation tensor. The small-scale noise is considered subsumed within the mesh and is not defined explicitly. This opens up various possibilities for turbulence modelling. The specification of the variance tensor a can be performed through an empirical local velocity fluctuation variance times a decorrelation time, or by physical models/approximations or using experimental measurements. The options explored in this study include physical approximation based models and empirical local variance based models as described below. Note that this derivation can be applied to any flow model. For instance, such a modelling has been successfully applied to derive stochastic large-scale representation of geophysical flows by (17; 33; 34).

A similar stochastic framework arising also from a decomposition of the Lagrangian velocity has been proposed in [START_REF] Holm | Variational principles for stochastic fluid dynamics[END_REF] and analysed in [START_REF] Crisan | Solution properties of a 3d stochastic Euler fluid equation[END_REF] and [START_REF] Cotter | Stochastic partial differential fluid equations as a diffusive limit of deterministic Lagrangian multi-time dynamics[END_REF]. This framework leads to enstrophy conservation whereas the formulation under location uncertainty conserves the kinetic energy of a transported scalar [START_REF] Resseguier | Geophysical flows under location uncertainty, part I: Random transport and general models[END_REF].

Physical approximation based models:

Smagorinsky's work on atmospheric flows and the corresponding model developement is considered to be the pioneering work on LES modelling [START_REF] Smagorinsky | General circulation experiments with the primitive equations[END_REF]. Based on Boussinesq's eddy viscosity hypothesis which postulates that the momentum transfer caused by turbulent eddies can be modelled by an eddy viscosity (ν t ) combined with Prandtl's mixing length hypothesis he developed a model (Smag) for characterising the SGS dissipation.

ν t = C||S||, (7)
τ = C||S||S, (8)
where τ stands for the SGS stress tensor, C is the Smagorinsky coefficient defined as (C s ∆) 2 , where ∆ is the LES filter width,

||S|| = 1 2 [ ij (∂ x i u j + ∂ x j u i ) 2 ] 1 2
is the Frobenius norm of the rate of strain tensor, and

S ij = 1 2 ( ∂ ūi ∂x j + ∂ ūj ∂x i ). (9)
Similar to Smagorinsky's eddy viscosity model, the variance tensor for the formulation under location uncertainty can also be specified using the strain rate tensor. Termed in the following as the Stochastic Smagorinsky model (StSm), it specifies the variance tensor similar to the eddy viscosity in the Classical Smagorinsky model:

a(x, t) = C||S||I 3 , (10)
where I 3 stands for 3 × 3 identity matrix and C is the Smagorinsky coefficient.

The equivalency between the two models can be obtained in the following case (as shown by ( 16)):

The SGS contribution (effective advection and SGS dissipation) for the StSm model is:

u j ∂ x j (∂ x j a kj ) + ij ∂ x i (a ij ∂ x j u k ) = u j ∂ x j (∂ x j ||S||δ kj ) + ij ∂ x i (||S||δ ij ∂ x j u k ), = u k ∆||S|| + ||S||∆u k + j ∂ x j ||S||∂ x j u k , (11)
and the SGS contribution for Smagorinsky model (∇ • τ ) is:

∇ • τ = j ∂ x j (||S||S), = j ∂ x j (||S||(∂ x j u k + ∂ x k u j )), = j ∂ x j ||S||∂ x j u k + ∂ x j ||S||∂ x k u j + ||S||∆u k . ( 12
)
An equivalency can be drawn between the two models by adding j ∂ x j ||S||∂ x k u ju k ∆||S|| to the StSm model. The additional term may also be written as:

∂ x k j ∂ x j (||S||)u j - j ∂ x j ∂ x k (||S||)u j -u k ∆||S||, ( 13
)
where the first term represents a velocity gradient which can be included within a modified pressure term as is employed for Smagorinsky model. The other two terms can be neglected for smooth enough strain rate tensor function. For smooth deformations both models are equivalent in terms of dissipation. It is important to note here that even if the effective advection is ignored in the StSm model, the two models still differ in a general case due to the first two terms in (13). Smagorinsky's pioneering work remains to date a popular model for LES, however it has certain associated drawbacks. The model assumes the existence of equilibrium between the kinetic energy flux across scale and the large scales of turbulence -this equilibrium is not established in many cases such as the leading edge of an airplane wing or turbulence with strong buoyancy. In addition, a form of arbitrariness and user dependency is introduced due to the presence of the Smagorinsky coefficient. This coefficient is assumed to be constant irrespective of position and time. Lilly [START_REF] Lilly | The representation of small scale turbulence in numerical simulation experiments[END_REF] suggests a constant coefficient value to be appropriate for the Smagorinsky model. However this was disproved by the works of [START_REF] Meyers | On the model coefficients for the standard and the variational multi-scale Smagorinsky model[END_REF], and ( 39) who shows that a constant value did not efficiently capture turbulence especially in boundary layers.

Numerous attempts were made to correct for the fixed constant such as damping functions [START_REF] Van Driest | The problem of aerodynamic heating[END_REF] or renormalisation group theory [START_REF] Yakhot | Renormalization group analysis of turbulence. I. Basic theory[END_REF]. Germano et al. [START_REF] Germano | A dynamic subgrid-scale eddy viscosity model[END_REF] provided a non ad-hoc manner of calculating the Smagorinsky coefficient varying with space and time using the Germano identity and an additional test filter t (termed as the Dynamic Smagorinsky (DSmag) model).

L ij = T ij -τ t ij , (14)
where τ stands for the SGS stress filtered by the test filter t, T is the filtered SGS stress calculated from the test filtered velocity field, and L stands for the resolved turbulent stress. The Smagorinsky coefficient can thus be calculated as:

C 2 s = < L ij M ij > < M ij M ij > , where (15)
M ij = -2∆ 2 (α 2 || St || St ij -(|| S|| Sij ) t (16)
and α stands for the ratio between the test filter and the LES filter. The dynamical update procedure removes the user dependancy aspect in the model, however it introduces unphysical values for the coefficient at certain instances. An averaging procedure along a homogenous direction is necessary to provide physical values for C s . However, most turbulent flows, foremost being wake flow around a cylinder, lack a homogenous direction for averaging. In such cases, defining the coefficient is difficult and needs adhoc measures such as local averaging, threshold limitation, and/or filtering methods to provide nominal values for C s . For the present study, we focus on the classical and dynamic variations of the Smagorinsky model -these model were used to study cylinder wake flow by ( 8), ( 42), [START_REF] Ouvrard | Classical and variational multiscale LES of the flow around a circular cylinder on unstructured grids[END_REF], among many others.

Local variance based models:

As the name states, the variance tensor can be calculated by an empirical covariance of the resolved velocity within a specified local neighbourhood. The neighbourhood can be spatially or temporally located giving rise to two formulations. A spatial neighbourhood based calculation (referred to as Stochastic Spatial Variance model (StSp)) is given as:

a(x, nδt) = 1 |Γ| -1 x i ∈η(x) (u(x i , nδt) -ū(x, nδt))(u(x i , nδt) -ū(x, nδt)) T C sp , (17)
where ū(x, nδt) stands for the empirical mean around the arbitrarily selected local neighbourhood defined by Γ. The constant C sp is defined as [START_REF] Harouna | Stochastic representation of the Reynolds transport theorem: revisiting large-scale modeling[END_REF]:

C sp = res η 5 3 ∆t, ( 18
)
where res is the resolved length scale, η is the Kolmogorov length scale and ∆t is the simulation time step. A similar local variance based model can be envisaged in the temporal framework; however, it has not been analysed in this paper due to memory limitations.

It is important to note that the prefix stochastic has been added to the MULU to differentiate the MULU version of the Smagorinsky model from its classical purely deterministic version. The model equations while derived using stochastic principles are applied in this work in a purely deterministic sense. The full stochastic formulation of MULU has been studied by [START_REF] Resseguier | Geophysical flows under location uncertainty, part I: Random transport and general models[END_REF].

Flow configuration and numerical methods

The flow was simulated using a parallelised flow solver, Incompact3d, developed by [START_REF] Laizet | High-order compact schemes for incompressible flows: A simple and efficient method with quasispectral accuracy[END_REF]. Incompact3d relies on a sixth order finite difference scheme (the discrete schemes are described in [START_REF] Lele | Compact finite difference schemes with spectral-like resolution[END_REF]) and the Immersed Boundary Method (IBM) (for more details on IBM refer to [START_REF] Gautier | A DNS study of jet control with microjets using an immersed boundary method[END_REF]) to emulate a body forcing. The main advantage of using IBM is the ability to represent the mesh in cartesian coordinates and the straightforward implementation of high-order finite difference schemes in this coordinate system. The IBM in Incompact3d has been applied effectively to cylinder wake flow by [START_REF] Parnaudeau | Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900[END_REF] and to other flows by [START_REF] Gautier | A DNS study of jet control with microjets using an immersed boundary method[END_REF], and (47) among others. A detailed explanation of the IBM as applied in Incompact3d, as well as its application to cylinder wake flow can also be found in [START_REF] Parnaudeau | Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900[END_REF]. It is important to note that this paper focuses on the accuracy of the sub-grid models within the code and not on the numerical methodology (IBM/numerical schemes) of the code itself.

The incompressibility condition is treated with a fractional step method based on the resolution of a Poisson equation in spectral space on a staggered pressure grid combined with IBM. While solving the Poisson equation for the stochastic formulation, the velocity bias was taken into account in order to satisfy the stochastic mass conservation constraints. It can be noted that although the solution of the Poisson equation in physical space is computationally heavy, the same when performed in Fourier space is cheap and easily implemented with Fast Fourier transforms. For more details on Incompact3d the authors refer you to [START_REF] Laizet | High-order compact schemes for incompressible flows: A simple and efficient method with quasispectral accuracy[END_REF] and [START_REF] Laizet | Incompact3d: A powerful tool to tackle turbulence problems with up to O(105) computational cores[END_REF].

The flow over the cylinder is simulated for a Re of 3900 on a domain measuring 20D × 20D × πD. The cylinder is placed in the centre of the lateral domain at 10D and at 5D from the domain inlet. For statistical purposes, the centre of the cylinder is assumed to be (0, 0). A coarse mesh resolution of 241 × 241 × 48 is used for the coarse LES (cLES). cLES discretisation has been termed as coarse as this resolution is ∼ 6.2% the resolution of the reference LES of (7) (henceforth referred to as LES -Parn). In terms of Kolmogorov units (η), the mesh size for the cLES is 41η × 7η -60η × 32η. The Kolmogorov length scale has been calculated based on the dissipation rate and viscosity, where the dissipation rate can be estimated as ∼ U 3 /L where U and L are the characteristic velocity scale and the integral length scale. A size range for y is used due to mesh stretching along the lateral (y) direction which provides a finer mesh in the middle. Despite the stretching, the minimum mesh size for the cLES is still larger than the mesh size of particle imagery velocimetry (PIV) reference measurements of (7) (henceforth referred to as PIV -Parn). For all simulations, inflow/outflow boundary condition is implemented along the streamwise (x) direction with free-slip and periodic boundary conditions along the lateral (y) and spanwise (z) directions respectively -the size of the spanwise domain has been fixed to πD as set by [START_REF] Beaudan | Numerical experiments on the flow past a circular cylinder at sub-critical Reynolds number[END_REF], which was also validated by [START_REF] Parnaudeau | Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900[END_REF] to be sufficient with periodic boundary conditions. The turbulence is initiated in the flow by introducing a white noise in the initial condition. Time advancement is performed using the third order Adam-Bashforth scheme. A fixed coefficient of 0.1 is used for the Smagorinsky models as suggested in literature (43) while a spatial neighbourhood of 7 × 7 × 7 is used for the Stochastic Spatial model. For the dynamic Smagorinsky model, despite the lack of clear homogenous direction, a spanwise averaging is employed. In addition, the constant is filtered and a threshold on negative and large positive coefficients is also applied to stabilise the model. Note that the positive threshold is mesh dependant and needs user-intervention to specify the limits.

The reference PIV [START_REF] Parnaudeau | Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900[END_REF] was performed with a cylinder of diameter 12 mm and 280 mm in length placed 3.5D from the entrance of the testing zone in a wind tunnel of length 100 cm and height 28 cm. Thin rectangular end plates placed 240 mm apart were used with a clearance of 20 mm between the plates and the wall. 2D2C measurements were carried out at a free stream velocity of 4.6 m s -1 (Re ∼ 3900) recording 5000 image pairs separated by 25 µs with a final interrogation window measuring 16 × 16 pixels with relatively weak noise. For more details about the experiment refer to [START_REF] Parnaudeau | Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900[END_REF].

The high resolution LES of ( 7) was performed on Incompact3d on the same domain measuring 20D × 20D × πD with 961 × 961 × 48 cartesian mesh points. The simulation was performed with the structure function model of ( 49) with a constant mesh size. LES -Parn is well resolved, however, there is a distinct statistical mismatch between LES -Parn and PIV-Parn especially along the centre-line (see figure 1a and figure 1b). Literature suggests that the wake behind the cylinder at a Re ∼ 3900 is highly volatile and different studies predict slightly varied profiles for the streamwise velocity along the centre-line. The averaging time period, the type of model, and the mesh type all affect the centre-line velocity profile. As can be seen in figure 1a and figure 1b, each reference data set predicts a different profile/magnitude for the streamwise velocity profiles. The DNS study of [START_REF] Ma | Dynamics and lowdimensionality of a turbulent near wake[END_REF] does not present the centreline velocity profiles. This provided the motivation for performing a DNS study at Re ∼ 3900 to accurately quantify the velocity profiles and to reduce the mismatch between the existing experimental and simulation datasets. The DNS was performed on the same domain with 1537 × 1025 × 96 cartesian mesh points using Incompact3d with stretching implemented in the lateral (y) direction.

From figure 1a we can see that the DNS and the PIV of Parnaudeau are the closest match among the data sets while significant deviation is seen in the other statistics. In the fluctuating streamwise velocity profiles, the only other data sets that exist are of (50) who performed Laser Doppler Velocimetry (LDV) experiments at Re = 3000 and Re = 5000. Among the remaining data-sets (LES of Parnaudeau, PIV of Parnaudeau, and current DNS) matching profiles are observed for the DNS and PIV despite a magnitude difference. These curves also match the profiles obtained by the experiments of Norberg (50) in shape, i.e. the similar magnitude, dual peak nature. The LES of Parnaudeau [START_REF] Parnaudeau | Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900[END_REF] is the only data-set to estimate an inflection point and hence is not considered further as a reference. The lower energy profile of the PIV may be attributed to the methods used for calculating the vector fields which employ a large-scale representation of the flow via interrogation win- 1 concisely depicts all the important parameters for the flow configuration as well as the reference datasets. Wake flow around a cylinder was simulated in the above enumerated configuration with the following SGS models: Classic Smagorinsky (Smag), Dynamic Smagorinsky (DSmag), Stochastic Smagorinsky (StSm), and Stochastic Spatial (StSp) variance. In accordance with the statistical comparison performed by [START_REF] Beaudan | Numerical experiments on the flow past a circular cylinder at sub-critical Reynolds number[END_REF], first and second order temporal statistics have been compared at 3 locations (x = 1.06D (top), x = 1.54D(middle), and x = 2.02D(bottom)) in the wake of the cylinder. All cLES statistics are computed (after an initial convergence period) over 90,000 time steps corresponding to 270 non-dimensional time or ∼ 54 vortex shedding cycles. All statistics are also averaged along the spanwise (z) direction. The model statistics are evaluated against the PIV experimental data of [START_REF] Parnaudeau | Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900[END_REF] and the DNS for which the data has been averaged over 400,000 time steps corresponding to 60 vortex sheddings. The work of [START_REF] Parnaudeau | Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900[END_REF] suggests that at least 52 vortex sheddings are needed for convergence which is satisfied for all the simulations. In addition, spanwise averaging of the statistics results in converged statistics comparable with the PIV ones. Both DNS and PIV statistics are provided for all statistical comparison, however, the DNS is used as the principal reference when an ambiguity exists between the two references.

Results

In this section, we present the model results, performance analysis and physical interpretations. Firstly, The cLES results are compared with the reference PIV and the DNS. The focus on centreline results for certain comparisons is to avoid redundancy and because these curves show maximum statistical deviation. This is followed by a characterisation and physical analysis of the velocity bias and SGS contributions for the MULU. The section is concluded with the computation costs of the different models.

Coarse LES

For cLES, the MULU have been compared with classic and dynamic version of the Smagorinksy model, DNS, and PIV -Parn. Figure 2 and figure 3 depict the mean streamwise and lateral velocity respectively plotted along the lateral (y) direction. In the mean streamwise velocity profile (see figure 2a), the velocity deficit behind the cylinder depicted via the U-shaped profile in the mean streamwise velocity is captured by all models. The expected downstream transition from U-shaped to a V-shaped profile is see for all the models -a delay in transition is observed for Smag model which biases the statistics at x = 1.54D and 2.02D. For the mean lateral component (see figure 3), all models display the anti-symmetric quality with respect to y = 0. Smag model shows maximum deviation from the reference DNS statistics in all observed profiles. All models but Smag capture the profile well while broadly StSp and DSmag models better capture the magnitude. As a general trend, Smag model can be seen to under-predict statistics while StSm model over-predicts. A better understanding of the model performance can be obtained through Figure 4 -6 which depict the second order statistics, i.e. the rms component of the streamwise (< u u >) and lateral (< v v >) velocity fluctuations and the cross-component (< u v >) fluctuations. The transitional state of the shear layer can be seen in the reference statistics by the two strong peaks at x = 1.06D in figure 4a. The magnitude of these peaks is in general underpredicted, however, a best estimate is given the MULU. DSmag and Smag models can be seen to under-predict these peaks at all x/D. This peak is eclipsed by a stronger peak further downstream due to the formation of the primary vortices (see figure 4b) which is captured by all the models.

The maxima at the centreline for figure 5 and the anti-symmetric structure for figure 6 are seen for all models. Significant mismatch is observed between the reference and the Smag/StSm models especially in figure 5a and 6a. In all second-order statistics, StSm model improves in estimation as we move further downstream. No such trend is seen for StSp or DSmag models while a constant under-prediction is seen for all Smag model statistics. This under-prediction could be due to the inherent over-dissipativeness of the Smagorinsky model which smooths the velocity field. This is corrected by DSmag/StSm models and in some instances over-corrected by the StSm model. A more detailed analysis of the two formulations under location uncertainty (StSm and StSp) is presented in sections 4.2.

The smoothing for each model is better observed in the 3D isocontours of vorticity modulus (Ω) plotted in figure 7. Plotted at non-dimensional Ω = 7, the isocontours provide an understanding of the dominant vortex structures within the flow. While large-scale vortex structures are observed in all flows, the small-scale structures and their spatial extent seen in the DNS are better represented by the MULU. The over-dissipativeness of the Smag model leads to smoothed isocontours with reduced spatial extent. The large-scale vortex structures behind the cylinder exhibit the spanwise periodicity observed by Williamson (53) for cylinder wake flow at low Re ∼ 270. Inferred to be due to mode B instability by Williamson, this spanwise periodicity was associated with the formation of small-scale streamwise vortex pairs. It is interesting to observe here the presence of similar periodicity at higher Re -this periodicity will be further studied at a later stage in this paper.

A stable shear layer associated with higher dissipation is observed in Smag model with the shear layer instabilities beginning further downstream than the MULU. An accurate shear layer comparison can be done by calculating the recirculation length (L r ) behind the cylinder. Also called the bubble length, it is the distance between the base of the cylinder and the point with null longitudinal mean velocity on the centreline of the wake flow. This parameter has been exclusively studied due to its strong dependence on external disturbances in experiments and numerical methods in simulations (54; 4). The effective capture of the recirculation length leads to the formation of U-shaped velocity profile in the near wake while the presence of external disturbances can lead to a V-shaped profile as obtained by the experiments of (9). Parnaudeau et al. [START_REF] Parnaudeau | Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900[END_REF] used this characteristic to effectively parameterise their simulations. The instantaneous contours can provide a qualitative outlook on the recirculation length based on shear layer breakdown and vortex formation. However, in order to quantify accurately the parameter, the mean and rms streamwise velocity fluctuation components were plotted in the streamwise (x) direction along the centreline (see figure 8a, and figure 8b). The recirculation length for each model is tabulated in table 2. StSp and DSmag models capture the size of the recirculation region with 0% error while the StSm model under estimates the length by 5.9% and the Smag model over estimates by 15.9%. The magnitude at the point of inflection is accurately captured by all the models (figure 8a).

For the rms centreline statistics of figure 8b, due to ambiguity between references, the DNS is chosen for comparison purposes. However, the similar magnitude, dual peak nature of the profile can be established through both the references. This dual peak nature of the model was also observed in the experiments of (50) who concluded that within experimental accuracy, the secondary peak was the slightly larger RMS peak as seen for the DNS. The presence of the secondary peak is attributed to the cross over of mode B for LES-Parn in figure 1b despite the simulation being within the transition regime.

The fluctuating centreline velocity profiles for the deterministic Smagorinsky models display an inflection point unlike the references. The MULU display a hint of the correct dual peak nature while under-predicting the magnitude matching with the PIV's large scale magnitude rather than the DNS. Although the Smag model has a second peak magnitude closer to the DNS, the position of this peak is shifted farther downstream. This combined with the inability of the model to capture the dual-peak nature speaks strongly against the validity of the Smag model statistics. Further analysis can be done by plotting 2D isocontours of the streamwise fluctuating velocity behind the cylinder, as shown in figure 9. The isocontours are averaged in time and along the spanwise direction. The profiles show a clear distinction between the classical models and the MULU in the vortex bubbles just behind the recirculation region. The vortex bubbles refer to the region in the wake where the initial fold-up of the vortices start to occur from the shear layers. The MULU match better with the DNS isocontours within this bubble as compared to the Smag or DSmag models. Along the centreline, MULU under-predict the magnitude, as depicted by the lower magnitude dual peaks in figure 8b. As we deviate from the centre-line, the match between the MULU and the DNS improves considerably. The mismatch of the isocontours in the vortex bubbles for the Smag and DSmag models with the DNS suggests that a higher magnitude for the centreline profile is not indicative of an accurate model. The dual peak nature of the streamwise velocity rms statistics show a strong dependance on the numerical model and parameters. This can be better understood via the constant definition within the StSp model formulation (refer to [START_REF] Harouna | Stochastic representation of the Reynolds transport theorem: revisiting large-scale modeling[END_REF]). The constant requires definition of the scale ( res ) of the simulation which is similar to ∆ used in classic Smagorinsky model, i.e. it defines the resolved length scale of the simulation. In the case of stretched mesh, the definition of this res can be tricky due to the lack of a fixed mesh size. It can be defined as a maximum (max(dx, dy, dz)) or a minimum (min(dx, dy, dz)) or an average (dx * dy * dz) (1/3) ). A larger value res is large, corresponding to a "PIV resolution", the centreline streamwise rms statistics display a dual peak nature with a larger initial peak similar to PIV reference. A smaller value for res , corresponding to a"higher resolution LES", shifts this dual peak into a small initial peak and a larger second peak similar to the DNS and of higher magnitude. The statistics shown above have been obtained with an res defined as (max(dx, dy, dz)) to emulate the coarseness within the model. Figure 10a and figure 10b show the power spectra of the streamwise and lateral velocity fluctuation calculated over time using probes at x/D = 3D behind the cylinder along the full spanwise domain. For the model power spectra, 135, 000 time-steps were considered corresponding to a nondimensional time of 405 which encapsulates ∼ 81 vortex shedding cycles. The hanning methodology is used to calculate the power spectra with an overlap of 50% considering 30 vortex shedding cycles in each sequence. The reference energy spectra (namely HWA) have been obtained from [START_REF] Parnaudeau | Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900[END_REF] while the DNS energy spectra has been calculated similar to the cLES. The process of spectra calculation for both reference and models are identical. All the values have been non-dimensionalised.

The fundamental harmonic frequency (f /f s = 1) and the second harmonic frequency are captured accurately by all models in the v-spectra. Let us recall that the cLES mesh is coarser than the PIV grid. Twice the vortex shedding frequency is captured by the peak in u-spectra at f /f s ∼ 2 as expected -twice the Strouhal frequency is observed due to symmetry condition at the centreline [START_REF] Ma | Dynamics and lowdimensionality of a turbulent near wake[END_REF]. The HWA measurement has an erroneous peak at f /f s ∼ 1 which was attributed to calibration issues and the cosine law by [START_REF] Parnaudeau | Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900[END_REF]. All models match with both the reference spectra. One can observe a clear inertial subrange for all models in line with the expected -5/3 slope. The models in the order of increasing energy at the small scales is DSmag <Smag = StSp <StSm. For the StSm model, an accumulation of energy is observed at the smaller scales in the u-spectra, unlike the StSp model. This suggests that the small-scale fluctuations seen in vorticity or velocity contours for the StSp model (i.e. in figure 7) are physical structures and not a numerical accumulation of energy at the smaller scales known to occur for LES.

The statistical comparisons show the accuracy and applicability of the MULU. The next sub-section focuses on the physical characterisation of the MULU -SGS dissipation, velocity bias and their contributions are studied in detail.

Velocity bias characterisation

The functioning of the MULU is through the small-scale velocity autocorrelation a. The effect of this parameter on the simulation is threefold: firstly, it contributes to a velocity bias/correction which is a unique feature of the MULU. Secondly, this velocity correction plays a vital part in the pressure calculation to maintain incompressibility. Finally, it contributes to the SGS dissipation similar to classical LES models -this signifies the dissipation occurring at small scales. This threefold feature of the MULU is explored in this section.

The contribution of the velocity bias can be characterised by simulating the MULU (StSm and StSp) with and without the velocity bias (denoted by Nad for no advection bias) and comparing the statistics (see figure 11a -11b). Only the centre-line statistics have been shown for this purpose as they display maximum statistical variation among the models and provide an appropriate medium for comparison. For the simulations without the velocity bias, the convective part in the NS equations remains purely a function of the large-scale velocity. In addition, the weak incompressibility constraint ( 5) is not enforced in the simulations with no velocity bias and the pressure is computed only on the basis of large-scale velocity. Similar to the Smagorinsky model, where the gradients of the trace of the stress tensor are considered subsumed within the pressure term, the divergence of the velocity bias is considered subsumed within the pressure term. The simulations parameters and flow configuration remain identical to cLES configuration.

The statistics show improvement in statistical correlation when the velocity bias is included in the simulation -all statistical profiles show improvement with inclusion of velocity bias but only the centre-line statistics have been shown to avoid redundancy. In the mean profile, the inclusion of velocity bias appears to correct the statistics for both models to match better with the reference. For the StSm model, there is a right shift in the statistics while the opposite is seen for the StSp model. The correction for the StSp model appears stronger than that for the StSm model. This is further supported by the fluctuation profile where without the velocity bias, the StSp model tends to the Smag model with an inflection point while the inclusion of a connection between the large-scale velocity advection and the small-scale variance results in the correct dual peak nature of the references. For the StSm model, figure 11b suggests a reduction in statistical correlation with the inclusion of the velocity bias -this is studied further through 2D isocontours.

Figure 12 plots 2D isocontours for the streamwise fluctuating profiles for the MULU. Once again an averaging is performed in time and along the spanwise direction. A clear distinction between the models with and without velocity bias is again difficult to observe. However, on closer inspection, within the vortex bubbles, we can see that including velocity bias improves the agreement with the DNS by reducing the bubble size for the StSm model and increasing it for the StSp model. The higher magnitude prediction along the centreline seen for the StSm -Nad model could be the result of an overall bias of the statistics and not due to an improvement in model performancethe presence of an inflection point in the profile further validates the model inaccuracy. This error is corrected in the model with velocity bias. This corrective nature of the bias is further analysed.

For the StSp model, the simulation without the bias has a larger recir-culation zone or is "over dissipative" and this is corrected by the bias. This result supports the findings of (55) whose structure model, when employed in physical space, applies a similar statistical averaging procedure of squarevelocity differences in a local neighbourhood. They found their model to also be over-dissipative in free-shear flows and did not work for wall flows as too much dissipation suppressed development of turbulence and had to be turned off in regions of low three-dimensionality. To achieve that [START_REF] Ducros | Large-eddy simulation of transition to turbulence in a boundary layer developing spatially over a flat plate[END_REF] proposed the filtered-structure-function model, which removes the large-scale fluctuations before computing the statistical average. They applied this model with success to large-eddy simulation and analysis of transition to turbulence in a boundary layer. For the StSp model, which also displays this over dissipation quality (without velocity bias), the correction appears to be implicitly done by the velocity bias. Such a velocity correction is consistent with the recent findings of ( 19) who provided physical interpretations of the local corrective advection due to the turbulence inhomegeneity in the pivotal region of the near wake shear layers where transition to turbulence takes place. The recirculation length for all cases is tabulated in table 3. Data are obtained from the centre-line velocity statistics shown in figure 11a. The tabulated values further exemplify the corrective nature of the velocity bias where we see an improved estimation of the recirculation length with the inclusion of the velocity bias. Also, a marginal improvement in statistical match similar to figure 11a is observed with the inclusion of the velocity bias for all lateral profiles (not shown here). It can be concluded that the inclusion of velocity bias provides, in general, an improvement to the model. The physical characteristics of the velocity bias (expressed henceforth as u * = 1 2 ∇ • a) are explored further. The bias u * , having the same units as velocity, can be seen as an extension or a correction to velocity. Extending this analogy, the divergence of u * is similar to "the divergence of a velocity field". This is the case in the MULU where to ensure incompressibility, the divergence free constraints (eq. ( 5)) are necessary. The stability and statistical accuracy of simulations were improved with a pressure field calculated using the modified velocity u, i.e. when the weak incompressibility constraint was enforced on the flow. This pressure field can be visualised as a true pressure field unlike the Smagorinsky model where the gradients of the trace of the stress tensor are absorbed in an effective pressure field.

Stretching the u * and velocity analogy, we can also interpret the curl of u * (∇ × u * ) as vorticity or more specifically as a vorticity bias. The curl of u * plays a role in the wake of the flow where it can be seen as a correction to the vorticity field. The divergence and curl of u * are features solely of the MULU and their characterisation defines the functioning of these models. Figure 13 depicts the mean isocontour of ∇ • (u * ) = 0.02 for the two MULU. This divergence function is included in the poisson equation for pressure calculation in order to enforce the weak incompressibility constraint. In the StSm model the contribution is strictly limited to within the shear layer while in the StSp model the spatial influence extends far into the downstream wake. The stark difference in the spatial range could be due to the lack of directional dissipation in the StSm model which is modelled on the classic Smagorinsky model. This modelling results in a constant diagonal auto correlation matrix, and the trace elements simplifying to a laplacian of a (∆a) for ∇ • (u * ). This formulation contains a no cross-correlation assumption (zero non-diagonal elements in the auto correlation matrix) as well as ignoring directional dissipation contribution (constant diagonal terms provides equal SGS dissipation in all three principle directions). These Smagorinsky like assumptions place a restriction on the form and magnitude of u * which are absent in the StSp model. The existence of cross-correlation terms in a for the StSp model results in a better defined and spatially well-extended structure for the divergence.

The importance of the cross-correlation terms are further amplified in the mean curl isocontour of u * (see figure 14) where once again a spatial limitation is observed for the StSm model. However, the more interesting observation is the presence of spanwise periodicity for the curl of u * observed in the StSp model. The curl parameter is analogous to vorticity and is in coherence with the birth of streamwise vortices seen in figure 7, a spanwise periodicity is observed with a wavelength λ ∼ 0.8. Figure 15 superimposes isocontour. While clear periodicity is not observed for the mean vorticity, alternate peaks and troughs can be seen which match with the peaks in the mean curl isocontour. The wavelength of this periodicity is comparable with the spanwise wavelength of approximately 1D of mode B instabilities observed by [START_REF] Williamson | Vortex dynamics in the cylinder wake[END_REF] for Re ∼ 270. The footprint of mode B instabilities is linked to secondary instabilities leading to streamwise vortices observed for Re ranging from 270 to ∼ 21000 [START_REF] Bays-Muchmore | On streamwise vortices in turbulent wakes of cylinders[END_REF]. These results demonstrate the ability of the spatial variance model to capture the essence of the auto-correlation tensor.

The regions of the flow affected by the auto-correlation term can be characterised by plotting the contours of SGS dissipation density ((∇u)a(∇u) T ) of the MULU averaged in time and along the spanwise direction. These have been compared with the dissipation densities for the Smag and DSmag models ((∇u)ν t (∇u) T ) (see figure 17). A 'reference' dissipation density has been obtained by filtering the DNS dissipation density at the cLES resolution (see figure 16e) and by plotting the difference. The StSp model density matches best with the DNS compared with all other models -a larger spatial extent and a better magnitude match for the dissipation density is observed. The high dissipation density observed just behind the recirculation zone is captured only by the StSp model while all Smag models under-predict the density in this region. The longer recirculation zone for Smag model can be observed in the density contours. A few important questions need to be addressed here: firstly, the Smag model is known to be over-dissipative, however, in the density contours, a lower magnitude is observed for this model. This is a case of cause and effect where the over-dissipative nature of the Smag model smooths the velocity field thus reducing the velocity gradients which inversely affects the value of the dissipation density. Secondly, in the statistical comparison only a marginal difference is observed especially between the DSmag and StSp models while in the dissipation density contours we observe considerable difference. This is because the statistical profiles are a result of contributions from the resolved scales and the sub-grid scales. The dissipation density contours of figure 17 represent only a contribution of the sub-grid scales, i.e. the scales of turbulence characterised by the model. Thus, larger differences are observed in this case due to focus on the scales of model activity. Finally, we observe in figure 9 that within the vortex bubbles behind the cylinder the MULU perform better than the Smag or DSmag models. For the StSp model, this improvement is associated with the higher magnitude seen within this region in the SGS dissipation density. For the StSm model, no such direct relation can be made with the SGS dissipation density. However, when we look at the resolved scale dissipation ((∇u)ν(∇u) T ) for the models (see figure 16), a higher density is observed in the vortex bubbles for this model. For the classical models high dissipation is observed mainly in the shear layers. As the kinematic viscosity (ν) for all models is the same, the density maps are indicative of the smoothness of the velocity gradient. For the classical models, we see a highly smoothed field while for the MULU, we see higher density in the wake. This difference could induce the isocontour mismatch seen in figure 9. These results are consistent with the findings of ( 19), who applied the MULU in the context of reduced order model and observed that MULU plays a significant role in the very near wake where important physical mechanisms take place. For the MULU, the SGS contributions can be split into velocity bias (u∇ T (-1 2 ∇ • a) and dissipation ( 1 However, it is important to note that while the computational cost for the Smagorinsky models stay fixed despite changes in model parameters, the cost for the StSp model strictly depends on the size of the local neighbourhood used. A smaller neighbourhood reduces the simulation cost but could lead to loss of accuracy and vice versa for a larger neighbourhood. The definition of an optimal local neighbourhood is one avenue of future research that could be promising. StSm model, which also provides comparable improvement on the classic Smag model, can be performed at 24% the cost of the DSmag model. Thus, the MULU provide a low cost (2/3rd reduction) alternative to the dynamic Smagorinsky model while improving the level of statistical accuracy.

2 ij ∂ x i (a ij ∂ x j u)).

Conclusion

In this study, cylinder wake flow in the transitional regime was simulated in a coarse mesh construct using the formulation under location uncertainty. The simulations were performed on a mesh 54 times coarser than the DNS study. The simulation resolution is of comparable size with PIV resolutionthis presents a useful tool for performing DA where in disparity between the two resolutions can lead to difficulties.

This study focused on the MULU whose formulation introduces a velocity bias term in addition to the SGS dissipation term. These models were compared with the classic and dynamic Smagorinsky model. The MULU were shown to perform well with a coarsened mesh -the statistical accuracy of the spatial variance based model was, in general, better than the other compared models. The spatial variance based model and DSmag model both captured accurately the volatile recirculation length. The 2D streamwise velocity isocontours of the MULU matched better with the DNS reference than the Smagorinsky models. Additionally, the physical characterisation of the MULU showed that the velocity bias improved the statistics -considerably in the case of the StSp model. The analogy of the velocity bias with velocity was explored further through divergence and curl functions. The spanwise periodicity observed at low Re in literature was observed at this higher Re with the StSp model through the curl of u * (analogous with vorticity) and through mean vorticity albeit noisily. The SGS contribution was compared with Smagorinsky models and the split of the velocity bias and dissipation was also delineated through isocontours.

The authors show that the performance of the MULU under a coarse mesh construct could provide the necessary computational cost reduction needed for performing LES of higher Re flows. The higher cost of the StSp model compared with Smagorinsky, is compensated by the improvement in accuracy obtained at coarse resolution. In addition, the StSp model performs marginally better than the currently established DSmag model at just 37% the cost of the DSmag model. This cost reduction could pave the way for different avenue of research such as sensitivity analyses, high Reynolds number flows, etc. Of particular interest is the possible expansion of Data Assimilation studies from the currently existing 2D assimilations [START_REF] Gronskis | Inflow and initial conditions for direct numerical simulation based on adjoint data assimilation[END_REF] or low Re 3D assimilations [START_REF] Robinson | Image assimilation techniques for Large Eddy Scale models : Application to 3d reconstruction[END_REF] to a more informative 3D assimilation at realistic Re making use of advanced experimental techniques such as tomo-PIV. Also and more importantly, the simplistic definition of the MULU facilitate an

Figure 1 :

 1 Figure 1: Mean streamwise velocity (a) and fluctuating streamwise velocity (b) in the streamwise direction along the centreline (y = 0) behind the cylinder for the reference data-sets. Legend: HWA -hot wire anemometry(7), K&M -B-spline simulations (case II) of Kravchenko and Moin(4), L&S -experiment of Lourenco and Shih (9), N -experiment of Norberg at Re = 3000 and 5000[START_REF] Norberg | LDV measurements in the near wake of a circular cylinder[END_REF], O&W -experiement of Ong and Wallace[START_REF] Ong | The velocity field of the turbulent very near wake of a circular cylinder[END_REF]

Figure 2 :

 2 Figure 2: Mean streamwise velocity at 1.06D (top), 1.54D (middle), and 2.02D (bottom) in the wake of the circular cylinder.

Figure 3 :

 3 Figure 3: Mean lateral velocity at 1.06D, 1.54D, and 2.02D in the wake of the circular cylinder.

Figure 4 :

 4 Figure 4: Streamwise rms velocity (u u ) fluctuations at 1.06D, 1.54D, and 2.02D in the wake of the circular cylinder.

Figure 5 :

 5 Figure 5: Lateral rms velocity (v v ) fluctuations at 1.06D, 1.54D, and 2.02D in the wake of the circular cylinder.

Figure 6 :

 6 Figure 6: Rms velocity fluctuations cross-component (u v ) at 1.06D, 1.54D, and 2.02D in the wake of the circular cylinder.

Figure 7 :

 7 Figure 7: 3D instantaneous vorticity iso-surface at Ω = 7.

Figure 8 :

 8 Figure 8: Mean (a) and Fluctuating (b) streamwise velocity profile in the streamwise direction along the centreline behind the cylinder.

Figure 9 :

 9 Figure 9: 2D isocontours of time averaged fluctuating streamwise velocity (u u ).

Figure 10 :

 10 Figure 10: Power spectra of streamwise (a) and lateral (b) velocity component at x/D = 3 behind the cylinder.

Figure 11 :

 11 Figure 11: Effect of velocity bias on centre-line mean (a) and fluctuating (b) streamwise velocity behind the cylinder.

Figure 12 :

 12 Figure 12: Effect of velocity bias on the 2D isocontour of time averaged fluctuating streamwise velocity (u u ).

  (a) Full scale view of the isocontour superposition with the outlined zoom area (b) Zoomed in view of the DNS mean vorticity isocontour (c) Zoomed in view of the isocontour superposition

Figure 15 :

 15 Figure 15: 3D isocontour superposition of the mean curl of u * (blue) for the StSp model at ∇ × (u * ) = 0.05 with the mean vorticity (yellow) for DNS at Ω = 3

Figure 16 :

 16 Figure 16: Sub-grid scale dissipation density in the wake of the cylinder. o stands for the original DNS dissipation and f stands for filtered (to cLES resolution) DNS dissipation.

Figure 18

 18 shows the contribution of the two via 3D isocontours (dissipation in yellow and velocity bias in red). The contribution of velocity bias is limited for the StSm model as expected while in the StSp model it plays a larger role. Velocity bias in StSp model is dominant in the near wake of the flow especially in and around the recirculation zone. It is important to outline that this is the

Figure 17 :

 17 Figure 17: Resolved scale dissipation density in the wake of the cylinder for each model.

Figure 18 :

 18 Figure18: 3d SGS contribution iso-surface along the primary flow direction (x) with dissipation iso-surface in yellow (at 0.002) and velocity bias in red (at 0.001)

  feature of the StSp model -the model captures accurately the statistics at only 0.37 the cost of performing the DSmag model.

Table 1 :

 1 Flow parameters. LES resolution[START_REF] Corpetti | Fluid experimental flow estimation based on an optical-flow scheme[END_REF]. The DNS, however, exhibits a profile similar to other references and a magnitude in between the two LDV experiments of Norberg. Considering the intermediate Reynolds number of the DNS compared to the Norberg experiments, this suggests good convergence and accuracy of the DNS statistics. Note that the cLES mesh is ∼ 0.46% the cost of the DNS. Table

		Re	nx × ny × nz	lx/D × ly /D × lz /D	∆x/D	∆y/D	∆z/D	U∆t/D
	cLES	3900	241×241×48	20×20×π	0.083	0.024-0.289	0.065	0.003
	DNS	3900	1537×1025×96	20×20×π	0.013	0.0056-0.068	0.033	0.00075
	PIV -Parn	3900	160×128×1	3.6×2.9×0.083	0.023	0.023	0.083	0.01
	LES -Parn	3900	961×961×48	20×20×π	0.021	0.021	0.065	0.003
	dows similar to a

Table 2 :

 2 Recirculation lengths for cLES.

	Model PIV -Parn DNS Smag DSmag StSm StSp
	L r /D	1.51	1.50 1.75	1.50	1.42 1.50
	of this parameter would signify a coarser mesh (i.e. a rough resolution) while
	a small value indicates a finer cut off scale or a finer mesh resolution. When

Table 3 :

 3 Recirculation lengths with and without velocity bias.

	Model PIV -Parn DNS StSm -Nad StSm StSp -Nad StSp
	L r /D	1.51	1.50	1.42	1.42	1.58	1.50


LES -SGS Models

Smagorinsky Variants :

Classic : ν t = (C s * ∆) 2 S,
(1)

Dynamic : C 2 s = - 1 2 L ij Mij M kl Mkl . (2)
where,

L ij = T ij -τij (3)
M ij = ∆2 | S| Sij -∆2 (| S| Sij ), (4)
WALE :

ν t = (C w ∆) 2 (ς d ij ς d ij ) 3/2 ( Sij Sij ) 5/2 + (ς d ij ς d ij ) 5/4
(5)

MULU

NS formulation as derived in [START_REF] Memin | Fluid flow dynamics under location uncertainty[END_REF] :

Mass conservation :

d t ρ t + ∇•(ρ w )dt + ∇ρ •σ dB t = 1 2 ∇• (a∇q)dt, (7)
w = w - 1 2 ∇ • a (8)
For an incompressible fluid :

∇• (σdB t ) = 0, ∇ • w = 0, (9)
Momentum conservation :

∂ t w +w ∇ T (w - 1 2 ∇ • a)- 1 2 ij ∂ x i (a ij ∂ x j w ) ρ = ρg -∇p+µ∆w . ( 10
)
Pranav CHANDRAMOULI 4D-Var with LES -Incompact3d Meeting 6 / 29

Modelling of a :

Stochastic Smagorinsky model (StSm) :

a(x, t) = C ||S||I 3 , (11)
Local variance based models (StSp / StTe) :

a(x, nδt) = 1 |Γ| -1 x i ∈η(x) (w (x i , nδt)-w (x, nδt))(w (x i , nδt)-w (x, nδt)) T C st , (
n x × n y × n z l x × l y ×

VDA -General Formulation

Objective : Estimate the unknown true state of interest x t (t, x)

Formulation :

∂ t x(t, x) + M(x(t, x)) = q(t, x), (13)
x(t 0 , x) = x b 0 + η(x) (14) Y(t, x) = H(x(t, x)) + (t, x) (15) q(t, x) -model error (Covariance matrix Q) η(x) -

VDA -General Formulation

Cost Function :

J(x 0 ) = 1 2 (x 0 -x b 0 ) T B -1 (x 0 -x b 0 )+ 1 2 t f t0 (H(x t ) -Y(t)) T R -1 (H(x t ) -Y(t)).

4DVar -Adjoint Method

Evolution of the state of interest (x 0 ) as a function of time

∂J ∂η = -λ(to) + B -1 (∂x(t 0 ) -∂x 0 ). ( 17
)
Le Dimet, Francois-Xavier, and Olivier Talagrand. "Variational algorithms for analysis and assimilation of meteorological observations : theoretical aspects." Tellus A : Dynamic Meteorology and Oceanography 38.2 (1986) : 97-110.

Pranav CHANDRAMOULI 4D-Var with LES -Incompact3d Meeting

Additional Control

Minimisation of J with respect to additional incremental control parameter δu :

δγ (i) = {δx (i) 0 , δu (i) } J(δγ (i) ) = 1 2 ||δx (i) 0 + x (i) 0 -x (b) 0 || 2 B -1 + 1 2 t f t 0 ||δu (i) + u (i) -u (b) || 2 Bc + ... (18)
constrained by :

∂ t δx (i) + ∂ x M(x (i) ) • ∂x (i) + ∂ u M(u (i) ) • ∂u (i) = 0 (19)
a(x, t) = C ||S||I 3 , (20)
Local variance based models (StSp / StTe) :

a(x, nδt) = 1 |Γ| -1 x i ∈η(x) (w (x i , nδt)-w (x, nδt))(w (x i , nδt)-w (x, nδt)) T C sp ,

τ

  Velocity fluctuation profiles for turbulent channel flow at Reτ = 395

  contour n x × n y × n z l x × l y × l z n y × n z l x /D × l y /D × l z /D ∆x/

Pranav CHANDRAMOULI4D-Var with LES -Incompact3d Meeting

Data Assimilation Introduction to Data Assimilation Code Formulation VDA with LES -Results Plan 1 Large Eddy Simulation SGS Models LES -Results 2 Data Assimilation Introduction to Data Assimilation Code Formulation VDA with LES -Results Data Assimilation Introduction to Data Assimilation Code Formulation VDA with LES -Results Data Assimilation Introduction to Data Assimilation Code Formulation VDA with LES -Results Data Assimilation -Types Sequential approach : Cuzol and Mémin (2009), Colburn et al. (2011), Kato and Obayashi (2013), Combes et al. (2015) Variational approach (VDA) : Papadakis and Mémin (2009), Suzuki et al. (2009), Heitz et al. (2010), Lemke and Sesterhenne (2013), Gronskis et al. (2013), Dovetta et al. (2014) Hybrid approach : Yang (2014) Pranav CHANDRAMOULI 4D-Var with LES -Incompact3d Meeting Data Assimilation Introduction to Data Assimilation Code Formulation VDA with LES -Results 4DVar -Incremental Optimisation Evolution of the cost function (J(x 0 )) as a function of time Courtier, P., J. N. Thépaut, and Anthony Hollingsworth. "A strategy for operational implementation of 4D-Var, using an incremental approach." Quarterly Journal of the Royal Meteorological Society 120.519 (1994) : 1367-1387. Pranav CHANDRAMOULI 4D-Var with LES -Incompact3d Meeting Data Assimilation Introduction to Data Assimilation Code Formulation VDA with LES -Results Code Formulation VDA with LES -Results 4DVar -Flow Chart 4DVar incremental Data Assimilation using Adjoint methodology. Pranav CHANDRAMOULI 4D-Var with LES -Incompact3d Meeting Data Assimilation Introduction to Data Assimilation Code Formulation VDA with LES -Results Data Assimilation Introduction to Data Assimilation Code Formulation VDA with LES -Results Data Assimilation Introduction to Data Assimilation Code Formulation VDA with LES -Results Data Assimilation Introduction to Data Assimilation Code Formulation VDA with LES -Results Data Assimilation Introduction to Data Assimilation Code Formulation VDA with LES -Results

Synthetic Assimilation at Re 3900

Optimisation Parameter -U(x, y, z, t 0 ) Control Parameters -U in (1, y, z, t) Cst (x, y , z)