Carlo Cintolesi
email: carlo.cintolesi@gmail.com

Etienne Mémin
email: etienne.memin@inria.fr

Stochastic modelling of turbulent flows for numerical simulations

Keywords: Stochastic models, Turbulence modelling, Numerical simulations, OpenFOAM

published or not. The documents may come

Stochastic modelling of turbulent flows for numerical simulations

Carlo Cintolesi, Etienne Mémin

Introduction

The reliable numerical simulation of turbulent flows is still nowadays a challenging issue, both in terms of mathematical modelling and of computational cost required. In the last decades, different techniques were developed to tackle this problem, the most fruitful for practical applications being the Reynoldsaveraged simulation and the Large-Eddy Simulation (LES) methodologies. Despite the continuous improvements with increasingF accuracy of the models, such methodologies are developed within a deterministic framework. Hence, they cannot completely represent the random nature exhibited by turbulent flows, that eventually requires the use of stochastic calculus. In the field of geophysical flows, probabilistic models are used to correct the effects of the coarse spatial discretisation. Similarly, the stochastic variables can be employed to account for the unresolved processes in the numerical reproduction of engineering and environmental flows.

The literature proposes different approaches on this topic. The stochastic Langevin equation is derived assuming that a fluid-particle velocity is perturbed by a Brownian motion, which is found to well described the dynamics of turbulent flows; see Pope [START_REF] Pope | Turbulent Flows[END_REF]. This equation was used in the framework of Probability Density Function (POF) methods to reproduce homogeneous isotropic turbulence, but also inhomogeneous and anisotropic turbulence by Pope [START_REF] Pope | A lagrangian two-time probability density function equation for inhomogeneous turbulent flows[END_REF] and by Durbin & Speziale [START_REF] Durbin | Realizability of second-moment closure via stochastic analysis[END_REF], respectively. Orszag [START_REF] Orszag | Analytical theories of turbulence[END_REF] and Leslie [START_REF] Leslie | Developments in the theory of turbulence[END_REF] introduced the Eddy-Damped Quasi-Normal Markovian (EDQNM) models; see the overview by Lesieur [START_REF] Lesieur | Turbulence in Fluids[END_REF]. The large-scale equations were closed in spectral space through a Gaussian closure. They were particularly suitable to study strong non-linearity in the small-scale turbulence. In the same framework, Chasnov [START_REF] Chasnov | Simulation of the kolmogorov inertial subrange using an improved subgrid model[END_REF] develops a forced-dissipative model, where the large-eddy Navier-Stokes equations were corrected by a stochastic force terms. This was a Gaussian forcing uncorrelated in time, homogeneous and isotropic in space.

Kraichnan [START_REF] Kraichnan | Dynamics of nonlinear stochastic systems[END_REF] exploits a different approach: the momentum equations are replaced by a set of equations with same mathematical properties, which are closed using a Gaussian stochastic model. This theory leads to valuable results in terms of mathematical properties (existence, singularities) and physical effects (turbulent diffusion, backscatter) analyses. Frederiksen [START_REF] Frederiksen | Subgrid modelling for geophysical flows[END_REF] shows that the same strategy can be used for a stochastic modelling of barotropic flows or in quasi-geostrophic approximation, that includes the interaction between topography and small-scale eddies. The randomness effects can be also explicitly introduced by means of ad hoc stochastic terms. Investigating the plane shear mixing layer, Leith [START_REF] Leith | Stochastic backscatter in a subgridscale model: Plane shear mixing layer[END_REF] improves the accuracy of LES with Smagorinsky model by introducing an explicit stochastic terms. On the theoretical side, Flandoli [START_REF] Flandoli | An Introduction to 3D Stochastic Fluid Dynamics[END_REF] studied fluid dynamic systems corrected with a random white noise force to reproduce the complex phenomena related to turbulence. These attempts have some limitations: the POF and EDQNM models required to work in the spectral space instead of the physical one; there is a certain degree of arbitrariness when explicit random terms are introduced (e.g. the random forcing should be multiplicative or additive); and overall the models can be hardly generalised for practical applications.

The methodology here presented aims to overcome these shortcomings. It develops from a different starting point: the fluid-particle trajectory in the Lagrangian framework is assumed to be a random process. It is expressed by a semimartingale, where the finite-variation part represents the smooth macroscopic velocity, while the martingale models the perturbations due to the turbulent motion. Consistently, an expression of the velocity is found and stochastic calculus is used to derived the stochastic equations of motions. In such a procedure, the use of the Itō-Wentzell formula is crucial to compute the time derivative, see Kunita [START_REF] Kunita | Stochastic Flows and Stochastic Differential Equations[END_REF]. A first work in this direction was that one of Brzeźniak [START_REF] Brzeźniak | Stochastic partial differential equations and turbulence[END_REF], subsequently extended by Mikulevicius and Rozovskii [27] and Flandoli [10]. Globally, these works focused on the mathematical properties of the stochastic equations. The work of Mémin [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF] follows a similar approach and developed the so called model under Location Uncertainty (LU), which is oriented to practical application in computational fluid dynamics. Recently, Holm [START_REF] Holm | Variational principles for stochastic fluid dynamics[END_REF] derived a similar set of equations using Lagrangian mechanics, which leads to additional terms, while Neves et al. [START_REF] Neves | Wellposedness for stochastic continuity equations with ladyzhenskaya-prodi-serrin condition[END_REF] studied theoretically a similar system of equations. The LU model was applied to different applications:

Resseguier et al. [START_REF] Resseguier | Geophysical flows under location uncertainty, Part I: random transport and general models[END_REF][START_REF] Resseguier | Geophysical flows under location uncertainty, Part II: quasi-geostrophy and efficient ensemble spreading[END_REF][START_REF] Resseguier | Geophysical flows under location uncertainty, Part III: sqg and frontal dynamics under strong turbulence conditions[END_REF] used it for geophysical flows simulations, where it exhibits a high accuracy in reproducing extreme events and provided new analysis tools. Chapron et al. [START_REF] Chapron | Largescale flows under location uncertainty: a consistent stochastic framework[END_REF] investigated the Lorentz-63 case and found that LU is able to explore the region of the deterministic attractor faster than the classical models. Resseguier et al. [START_REF] Resseguier | Stochastic modelling and diffusion modes for proper orthogonal decomposition models and small-scale flow analysis[END_REF] employed it in conjunction with the proper orthogonal decomposition technique for the numerical simulation of a flow past a circular cylinder at Re = 3900.

Although this is a promising methodology, the inherent mathematical complexity of stochastic partial differential equations poses some difficulties: the resolution of stochastic partial equations is not straightforward and can considerably increase the simulation time. For these reasons, Mémin [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF] also introduces a simplified model, where the resolution of stochastic equations is avoided by modelling the effects of the random velocity term by physical assumption. This give rise to the so called pseudo-stochastic simulation (PSS) methodology: the flow dynamics is described by classical partial differential equations, which includes additional terms provided by the stochastic modeling. The PSS was adopted by Harouna and Mémin [START_REF] Harouna | Stochastic representation of the reynolds transport theorem: Revisiting large-scale modeling[END_REF] to investigate the Green-Taylor vortex flow applying several models for the stochastic contribution. Chandramouli et al. [START_REF] Chandramouli | Coarse large-eddy simulations in a transitional wake flow with flow models under location uncertainty[END_REF] employed it to simulate the transitional wake flow with coarse mesh resolution, proving that it generates a more accurate outcomes with respect to classical LES.

Notwithstanding the above mentioned studies, a pointwise analysis of the pseudo-stochastic model is lacking. The aims of the present work is to study in details the characteristics of the LU and the PSS model, both theoretically and numerically, establishing a parallelism with the classical LES methodology. First, a theoretical analysis of the PSS equations is reported; second, a simplified closure model is adopted to perform numerical simulations on the plane channel flow at Re τ = 590. The simulation outcomes are discussed in light of the previous theoretical analysis and the peculiarity of the PSS are highlighted. The main novelty of this work is to propose a detailed and systematic comparison between PSS and LES approach, pointing out the physical meaning of the extra term arising from the stochastic derivations (supported by simulations). Moreover, after few years from its first formulation, an alternative mathematically derivation of the LU and PSS model is proposed in the appendix. Efforts have been made to simplify and give a linear structure to the procedure, highlighting the key hypotheses.

The paper is organized as follows: section 2 presents the pseudo-stochastic model and the relative turbulent kinetic budget; section 3 reports a physical interpretation of the model and make a comparison with the LES methodology; section 4 describe the closure model for PSS; section 5 discusses the numerical simulation results; section 6 gives some final remarks. In appendix A an alternative and detailed derivation of the stochastic model for turbulent flows is presented.

Pseudo-Stochastic Model

In this section, the stochastic formalism and the pseudo-stochastic equations are reported.

Stochastic formalism

The particle trajectory in a turbulent regime is not completely known because it is subject to some random (turbulent) effects. Consequently, the fluidparticle displacement is described by the stochastic differential equations of the type:

dX i t (x 0 ) = w i (X t , t)dt + dη i t (X t ), (1)
where the index i = 1, 2, 3 indicates respectively the x,y,z-component in space (they are placed at top or bottom indifferently); X i t (x 0 ) is the trajectory followed by a fluid-particle initially located in x 0 ; w i is a differentiable function of bounded variation (i.e. equivalent to a deterministic function) that corresponds to the resolved flow velocity; η i t = t 0 dη i t is a martingale that accounts for the stochastic contributions to the motion. The Einstein summation convention over repeated indexes is adopted. The stochastic contribution is constructed as a combination of a cylindrical Wiener processes B k t (x) not differentiable in time, and a differentiable diffusion tensor σ ik which acts as an integral kernel:

dη i t (x) = Ω σ ik (x, y, t)dB k t (y) dy. (2)
Notice that the stochastic processes η i t are uncorrelated in time and correlated in space by means of the diffusion tensor.

The expression of the velocity field U i in Eulerian coordinate x is derived from equation (1); it reads:

U i (x, t) = w i (x, t) + ηi t (x), (3)
where the second term on the right-hand side expresses the stochastic velocity defined by formula [START_REF] Leslie | Developments in the theory of turbulence[END_REF]. From a physical point of view, w i is the velocity expected value and ηi t (x) represents a noise: a generalised stochastic process that has to be defined in the space of temperate distribution, see Øksendal [START_REF] Øksendal | Stochastic Differential Equations[END_REF].

The quadratic variation of the diffusion tensor is a quantity of particular interest; it represents the time-variation of the spatial variance of the stochastic increments along time. It is named variance tensor and is defined as:

a ij (x, t) = Ω σ ik (x, y, t)σ jk (x, y, t) dy. (4)
As a function, it is assumed to have all the regularity (differentiable and integrable in time and space) required by computation; as a tensor, it is a pointwise symmetric and semi-positive definite matrix.

Pseudo-stochastic equations of motion

The stochastic fluid dynamics equations for a Newtonian incompressible fluid are derived in appendix A. The final system (63) is composed by one set of stochastic equations and one of pure deterministic ones. The former allows to find an expression for the variance tensor a ij , which is required for the resolution of the latter. Together, they provide a close system of equations that composes the LU model. Let us not that full stochastic model can be obtained by relaxing the assumption of bounded variation for the resolved velocity (see [START_REF] Resseguier | Geophysical flows under location uncertainty, Part I: random transport and general models[END_REF]).

In order to simplify the model by avoiding the resolution of stochastic partial differential equations, the variance tensor a ij is not computed but modelled through physical assumptions. This choice gives rise to a hybrid model where the stochastic contribution on the governing equations is modelled by a deterministic function, and, overall, no stochastic equations have to be resolved.

Such model leads to pseudo-stochastic simulation approach. The PSS momentum and continuity equations for incompressible flows reads, respectively:

       ∂w i ∂t + w * j ∂w i ∂x j = - ∂p ∂x i + ν ∂ 2 w i ∂x j ∂x j + 1 2 ∂ ∂x j a jk ∂w i ∂x k ∂w * i ∂x i = 0. ( 5
)
where ν is the molecular viscosity, the modified pressure p = p h + ν 3 ∂w ∂x is the sum of the hydrostatic pressure and the divergence of the velocity field (which is not solenoidal), and the effective advection velocity w * i reads:

w * i = w i - 1 2 ∂ ∂x k a ik . (6)
The terms depending on a ij account for the effects of the Stochastic Unresolved Scales (SUS) of motion, since the variance tensor is a measure of the intensity and the anisotropy of turbulent random velocities.

Notice that system (5) reduces to the classical Navier-Stokes equations when the a ij is the zero matrix, i.e. when the stochastic contributions disappear.

Resolved kinetic energy budget

The turbulent kinetic energy (TKE) budget of the resolved scales of motion is presented. The resolved velocity is decomposed in a mean and a fluctuating part, respectively:

w i = W i + w i , (7)
where the capitol letter indicates the averaged field, W i = w i . Variance tensor and pressure are decomposed in a similar way: a ij = A ij + a ij and p = P + p . The (resolved) turbulent kinetic energy κ = w i w i /2 budget reads:

∂ κ ∂t + W j - ∂ ∂x k A jk 2 ∂ κ ∂x j + w j - ∂ ∂x k a jk 2 ∂κ ∂x j advection = = ∂ ∂x j -p w j + νδ jk + A jk 2 ∂ κ ∂x j + a jk 2 ∂κ ∂x j + a jk w i 2 ∂W i ∂x k transport + p 2 ∂ 2 a jk ∂x j ∂x k turb. compress.
νδ jk + A jk 2

∂w i ∂x j ∂w i ∂x k - a jk 2 ∂w i ∂x j ∂w i ∂x k dissipation -w j - ∂ ∂x k a jk 2 w i ∂W i ∂x j production - a jk 2 ∂w i ∂x j ∂W i ∂x k loss to SU S (8)
The TKE terms are interpreted in light of the classical budget analysis, e.g. see Kundu and Cohen [START_REF] Kundu | Fluid Mechanics[END_REF]. On the left-hand side, the second and third terms represent the TKE advection by mean and SUS effective advection velocity.

On the right-hand side:

first four terms: transport by pressure, molecular viscosity and turbulent stresses;

fifth term: turbulent compression/expansion due to SUS;

sixth and seventh terms: dissipation by molecular viscosity (it can be proven that A ij is positive defined), resolved turbulence and SUS motions;

eight term: shear production, this term appears in the mean kinetic budget (not shown here) with opposite sign;

last term: loss due to SUS also present in the mean kinetic energy budget.

The pseudo-stochastic TKE budget reduces to the classical one if the stochastic contribution is negligible a ij 0. It is worth to notice that the production term includes the contribution of the fluctuations of turbulent advection velocity, while the variance tensor plays a role of a turbulent viscosity dissipation tensor.

Analysis of pseudo-stochastic model

The expression of fluid-particle displacement (1) states that a particle trajectory is driven by two actors: a differentiable drift velocity and a Brownian process highly fluctuating in time. In the framework of PSS, the drift velocity w i that can be interpreted as the resolved velocity field, while the random field assembles the residual motion that are fast oscillating stochastic components, possibly anisotropic and non-homogeneous in space.

Physical interpretation

Recalling the decomposition of the velocity gradient in symmetric and antisymmetric parts, respectively called the strain-rate tensor and the rotation-rate tensor:

∂w i ∂x j = 1 2
∂w i ∂x j + ∂w j ∂x i + 1 2 ∂w i ∂x j - ∂w j ∂x i = S ij + Ω ij , (9)
the pseudo-stochastic Navier-Stokes equation ( 5) and continuity equation ( 6) are rearranged as, respectively:

∂w i ∂t + w j - 1 2
∂a jk ∂x k ef f ective advection ∂w i ∂x j = - ∂ ∂x i p h + ν 3 ∂ 2 a sk ∂x k ∂x s modif ied pressure +2ν ∂S ij ∂x j + 1 2 ∂ ∂x s (a sk S ki ) - 1 2 ∂ ∂x s (a sk Ω ki ) dif f usion due to SU S , (10)
and

∂w i ∂x i = 1 2 ∂ 2 a jk ∂x j ∂x k turb. compr. . (11)
The terms that depend on variance tensor account for the influence of the SUS on the resolved scales. A physical interpretation of such terms is proposed:

Effective advection: the advection velocity is corrected by an inhomogeneous turbulence contribution. It corresponds to a velocity induced by the unresolved turbulent motions, that can be linked to the turbophoresis phenomenon detectable in geophysical flows; i.e. the tendency of fluid-particle to migrate in the direction of less energetic turbulence (see also [START_REF] Resseguier | Geophysical flows under location uncertainty, Part I: random transport and general models[END_REF]).

Modified pressure: the non-solenoidal velocity field leads to the presence of an isotropic turbulent factor, that has the dimension of a pressure:

p t = ν 3 ∂ 2 a sk
∂x k ∂xs . This term does not contribute to the flow and it is included in the pressure gradient in the same manner as the isotropic residual stress in the Smagorinsky model, see [START_REF] Pope | Turbulent Flows[END_REF].

Diffusion due to SUS: they account for the turbulent diffusion; the variance tensor plays the role of a diffusion tensor similar to a generalised eddy-viscosity matrix. Both the deformation rate and rotation-rate contribute to diffusion, unlike to the classical eddy-viscosity model in which fluid rotation-rate is assumed to be irrelevant in turbulent modelling (see also following section 3.2).

Turbulent compressibility: the continuity equation [START_REF] Frederiksen | Subgrid modelling for geophysical flows[END_REF] suggests that the flow is turbulent-compressible; i.e. the unresolved turbulence induces a local fluid compression or expansion.

The variance tensor is the key parameter of the pseudo-stochastic model. It has the physical dimension of a dynamic viscosity [m 2 /s] and carries information on the intensity and the anisotropy of the SUS. As already mentioned, a ij can be interpreted as a generalised eddy-viscosity parameter. Implicitly, this leads to the hypothesis that the SUS influences the resolved flow as an alteration of fluid viscosity, that is an empirical consideration largely accepted.

The divergence of the variance tensor is hereafter named turbulent advection velocity:

u ta,i = - 1 2 ∂a ij ∂x j , (12)
while the divergence of the turbulent advection velocity measured the turbulent compressibility:

Φ tc = 1 2 ∂ 2 a ij ∂x i ∂x j , (13)
and it is directly proportional to the isotropic turbulent factor p t appearing in the modified pressure. The numerical simulations reported later allow to gain additional insights regarding these two quantities, we refer to section 5.3 for the numerical analysis.

Comparison with LES eddy-viscosity models

The LES methodology consists in applying a spatial filter to velocity field, and then directly resolve the filtered velocity and model the sub-filter velocities.

See Sagaut [START_REF] Sagaut | Large eddy simulation for incompressible flows[END_REF] and Piomelli [START_REF] Piomelli | Large-eddy and direct simulation of turbulent flows[END_REF] for an extended introduction on this subject.

Practically, the computational grid acts as an implicit spatial filter on the governing equations, which generates an extra term τ ij in the classical Navier-Stokes equations:

       ∂u i ∂t + u j ∂u i ∂x j = - ∂p ∂x i + ν ∂ 2 u i ∂x j ∂x j - ∂τ ij ∂x j , ∂u i ∂x i = 0, (14)
where the sub-grid scale (SGS) tensor is τ ij = u i u ju i u j , and the straight over-bar denotes the spatial filter associated to the local cell width, computed as ∆ = (∆x∆y∆z) 1/3 . Adopting the eddy-viscosity assumption, the anisotropic part of such tensor reads:

τ R ij = τ ij - τ kk 3 δ ij = -2ν sgs S ij , (15)
where ν sgs is the SGS viscosity parameter, which has to be specified by additional models (e.g. Smagorinsky model, Spalart-Allmaras, kω, k -). Equation (15) implies that: (a) the anisotropic Reynolds stress tensor is aligned with the mean strain-rate tensor; (b) the two are directly proportional through a single parameter, equal for all the six independent components of τ R ij .

The pseudo-stochastic model is equivalent to a constant eddy-viscosity model if the variance tensor is expressed by a ij = 2ν sus δ ij where the SUS viscosity ν sus is constant. In this sense, the pseudo-stochastic model can be considered a generalisation of the eddy-viscosity model. The theoretical advantages of the former to the latter are pointed out:

1. The PSS does not rely on hypothesis (a). The effects of unresolved scales of motion are given by a ij , without imposing any constrains on the directions along with the SUS acts on the resolved flow. The eddy-viscosity models are quite reasonable for simple shear flows and it is largely applied in computational fluid dynamics. However, most of their shortcomings derive from the fact that hypotheses (a) and (b) are not generally satisfied; see Pope [START_REF] Pope | Turbulent Flows[END_REF]. Efforts have been made to develop alternative models where the principal axis of τ R ij are not forced to be aligned with those of the mean strain tensor (e.g. the Reynolds-stress models), or where equation ( 15) is substituted by a non-linear viscosity models, in which the rotation strain-rate comes into play, see for example Bauer et al. [START_REF] Bauer | Accuracy and robustness of nonlinear eddy viscosity models[END_REF]. In geophysical flow simulations, the strong grid anisotropy between horizontal and vertical directions is successfully handled using a directional eddy-viscosity, see Roman

and Armenio [START_REF] Roman | Large eddy simulation of mixing in coastal areas[END_REF].

It is worth mentioning that the eddy-viscosity parameter a ij comes directly from the basic assumption of velocity decomposition in a smooth and a fast oscillating components (3), whereas it is introduced in LES equations through an ad hoc physical assumption.

Variance tensor model

In the LES framework, a popular model for ν sgs in LES methodology is the the Smagorinsky model, first proposed by Smagorinsky [START_REF] Smagorinsky | General circulation experiments with the primitive equations: I. the basic experiment[END_REF] for simulation of environmental flows (see also Deardorff [START_REF] Deardorff | A numerical study of three-dimensional turbulent channel flow at large reynolds numbers[END_REF]). It is derived under the hypothesis of local equilibrium between production and dissipation of turbulent kinetic energy, and reads:

ν sgs = c 2 s ∆ 2 |S|, ( 16
)
where |S| is the norm of the strain-rate tensor. The parameter c 2 s is set constant and can be evaluated from experiments, direct numerical simulations or analytical considerations, e.g. see Lilly [START_REF] Lilly | The representation of small-scale tubulence in numerical simulation experiments[END_REF].

In order to perform a close comparison with the LES methodology, the variance tensor is modelled by a simple model analogous to the Smagorinsky model:

a ij = c m ∆ 2 |S|δ ij ( 17
)
where ∆ is the cell grid width and c m is a model parameter. Hence, the variance tensor reduces to a diagonal matrix with equal elements because turbulence is assumed isotropic and homogeneous in all directions.

The relation with the classical Smagorinsky model is now highlighted. In LES, having applied the Smagorinsky model, the anisotropic Reynolds stress tensor reads:

- ∂τ R ij ∂x i = ∂ ∂x i 2c 2 s ∆ 2 |S|S ij = S ij ∂C s |S| ∂x i + C s 2 |S| ∂ 2 w j ∂x i ∂x i , (18)
where C s = 2c 2 s ∆ 2 denotes an auxiliary variable, c 2 s is the Smagorinsky parameter and the velocity is divergence-free. In the PSS, the total turbulent model can be expressed by a single term, that gathers the dissipative and turbulent advective contributions. Applying formula [START_REF] Kundu | Fluid Mechanics[END_REF] with c m = 2c 2 s and defining

C m = c m ∆ 2 , such a term becomes: 1 2 ∂ 2 a sk w i ∂x s ∂x k = S ij ∂C m |S| ∂x i + C m 2 |S| ∂ 2 w j ∂x i ∂x i eddy-viscosity terms + Ω ij ∂C m |S| ∂x j rotational term + w i 2 ∂ 2 C m |S| ∂x j ∂x j strain-rate dif f. , (19)
where the first two terms on the right-hand side have (formally) the same expression as [START_REF] Kunita | Stochastic Flows and Stochastic Differential Equations[END_REF], while the third and fourth term are additional contributions.

The PSS with isotropic constant model reduces to the LES Smagorinsky model under two approximations:

1. the rotation-rate does not contribute to turbulence effects on the mean flow, thus it is neglected;

2. the norm of strain-rate tensor is almost harmonic (Laplacian is close to zero), which makes the fourth term negligible.

Notice that the latter hypothesis implies that the continuity equation ( 6 The stretching is symmetric with respect to the channel centre plane y = δ, and it is obtained with a double-side stretching function based on hyperbolic tangent:

y(ξ) = 1 2 1 + tanh(λ(ξ -1/2)) tanh(λ/2) , ( 20
)
where ξ is the vertical coordinate of uniform point distribution and the stretch- & Issa [START_REF] Oliveira | An improved piso algorithm for the computation of bouyancy driven flows[END_REF].

Variables are discretised in space with a second-order central difference scheme, while time integration is performed using an implicit Euler backward scheme. Such a scheme employs the variables at the previous two time steps, leading to a second order accuracy. Globally, the numerical solvers are secondorder accurate in time and space. The time advancement fulfils the Courant-Friedrichs-Lewy condition Co < 0.5. The Courant number is computed as Co = ∆t|v|/δx, where: ∆t is the time step, |v| is the velocity magnitude through the cell, δx is the cell length. The model constants are chosen to be c 2 s = c m /2 = 0.004225, and for PSS v = w * while for LES w * = u.

Results discussion

The simulations are run till the statistical steady state is reached, then they are re-run for an additional period of 12t 0 where the statistics are collected.

The quantities are averaged in time and in space along span-wise and streamwise directions, and exploiting the domain symmetry in vertical direction. The angular brackets ψ denote the average in time and wall-parallel directions for a generic variable ψ.

First and second order statistics

The first and second order statistics of the velocity filed are analysed.

Effects of the extra terms in PSS

The LES constant Smagorinsky model and the PSS constant isotropic model lead to similar governing equations, but the latter has some additional terms not present in the former: the eddy-viscosity terms (EV), the rotational term (RT) and the strain-rate diffusion (SD) defined in equation [START_REF] Gall | Brownian Motion, Martingales, and Stochastic Calculus[END_REF]. The influence of such terms is checked.

Figure 2 shows the x-component of the above-mentioned terms (averaged)

versus the wall coordinate. They are made non-dimensional by u 3 τ /ν. In LES, the term EV accounts for all sub-grid scale effects and represents a negative turbulent diffusion near the wall. In the PSS constant isotropic model, two other terms come into play: SD is negative in the region y + < 10, while it

shows positive value at y + > 10; RT exhibits a positive contribution against the negative one of EV. The three terms become negligible in the log-law region; hence, the SUS model acts mainly at the near-wall region. The point y + = 10, located in the buffer layer, is of particular interest: approximately at this height, EV and RT reach the minimum and maximum (respectively), while SD changes sign. Globally, the RT and SD terms reduce the negative contribution of EV to the velocity equations in the buffer region, eventually producing a positive turbulent diffusion. negative values, with a climax at y + ∼ = 10. Quantitatively, the turbulent advection is not strong enough to produce remarkable results on the mean flow; however, it generates a weak vertical velocity w y directed from the center to the wall of the channel (not reported). Hence, u ta is qualified as a weak turbophoresis velocity: it advects the flow from the buffer region to the log-law region, i.e. in the direction of decreasing turbulence level (estimated by the velocity RMS intensity). The turbulent compressibility Φ tc assumes negative values in the viscous sub-layer and positive values in the buffer layer. Elsewhere, it is practically zero. In light of equation [START_REF] Frederiksen | Subgrid modelling for geophysical flows[END_REF], this behaviour is related to the presence of a turbulent fluid compression and expansion, respectively.

Turbulent advection and compressibility

Additional insight on this phenomenon is gained visualising the Φ tc instantaneous values. same stream-wise extension and span-wise period. Also, the negative spots are located at the same height at which steaks are triggered. Therefore, these two structures appears to be related.

Resolved turbulent kinetic budget

The pseudo-stochastic TKE budget ( 8) is finally scrutinised for PSS and LES simulations Figure 5 shows selected terms of the TKE budget. Production and dissipation profiles are similar for PSS and LES, but the former appears to be more effective in energy dissipation in near-wall region and has a higher production of TKE in the range 5 < y + < 20. In the PSS, the turbulent compression term is almost zero and does not contribute to the budget; while the loss due to SUS presents slightly negative values mainly localised in the viscous layer.

Hence, it contributes to global energy dissipation.

Conclusions

The pseudo-stochastic simulation (PSS) methodology introduced by Mémin [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF] is analysed theoretically and numerically, through a direct comparison with The PSS of turbulent channel flow at Re τ = 590 is performed, together with the LES with constant Smagorinsky sub-grid scale model. For a better comparison, a closure model analogous to the Smagorinsky one is used for the PSS. However, it is shown that this last does not rely on the eddy-viscosity hypotheses, hence it is not affected by its shortcomings. The PSS does not show improvement in first and second order statistics, possibly because of the simple expression of a ij , but reproduces additional features: a weak turbophoresis is detected in the buffer region, while a turbulent compression and expansion is identified in the viscous and buffer layer (respectively). This quantity appears to be related to the streaks, turbulent structures appearing near the wall region.

Finally, the pseudo-stochastic model is a generalisation of the LES eddyviscosity model and describes a richer physics. Overall, it represents a promising approach for simulation of turbulent flows: the mathematical analysis here reported gives a clear physical interpretation of the model, supported by numerical results.

A . Formal derivation of stochastic model

The mathematical conditions under which this derivation is consistent are reported in [START_REF] Flandoli | An Introduction to 3D Stochastic Fluid Dynamics[END_REF][START_REF] Mikulevicius | On Equations of Stochastic Fluid Mechanics[END_REF]27]. An introduction to the mathematical framework in which the present model is developed can be found in Øksendal [START_REF] Øksendal | Stochastic Differential Equations[END_REF] and Kunita [START_REF] Kunita | Stochastic Flows and Stochastic Differential Equations[END_REF].

A.1 Trajectory and stochastic velocity definitions

As already mentioned, expression (1) has to be understood in an integral sense:

X i t = X i 0 + t 0 w i (Xs, s)ds + t 0 dη i s (Xs), (21)
where the Itō stochastic integral is used to integrate the random process. The process

X i t (x 0 )
is a semimartingale defined for each spatial point x 0 ∈ Ω and time t ∈ T ⊆ R + in an appropriate probability space.

The stochastic velocity in equation ( 3) is a symbolic expression that is defined as a weak derivative of the random displacement:

h(t) ηi t (x)dt = h (t)η i t (x)dt, (22)
for each h test function; see also [START_REF] Øksendal | Stochastic Differential Equations[END_REF].

A.2 The stochastic Reynolds transport theorem

Being the velocity field a stochastic process, the governing equations of fluid dynamics cannot be recovered using deterministic calculus, ref. [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF]27]. In this concern, the key point is to give an expression of the Reynolds transport theorem (RTT) for stochastic quantities.

Subsequently, the stochastic Navier-Stokes equations are found imposing conservation of mass and momentum.

Theorem 1 (Stochastic RTT) Let us consider a physical quantity q(x, t) within a material volume V (t) ⊂ R 3 , transported by a stochastic flow of the form (1) and such that it can be written as a semimartingale of the type:

q(x, t) = q(x, 0) + t 0 g(x, s)ds + t 0 Ω f k (x, y, s)dB k s (y)dyds, (23)
where g, f are processes of bounded-variation and the Itō integral are employed. If the following properties holds:

1. symmetric diffusion tensor:

σ ij = σ ji , 2. solenoidal diffusion tensor: ∂ ∂x i σ ij (x, y, t) = 0 for all j,
3. conserved quantity: dq(Xt, t) = 0, then the stochastic RTT has an explicit differential form that reads:

d V (t) q(x, t)dx = V (t) ∂tq + ∂(qw i ) ∂x i dt - 1 2 ∂ 2 (qa ij ) ∂x i ∂x j dt + ∂q ∂x i dη i t dx, ( 24
)
where ∂t is the differential with respect to the second variable, and d denotes the total time increment at a fixed spatial point.

Assumption (1) greatly simplifies the computation and can be justified a posteriori: it implies that the variance tensor a ij is symmetric, a desirable properties in light of its physical interpretation (see section 3.1). Therefore, this assumption is considered reasonable in the fluid dynamics context; however it is not mandatory, see [START_REF] Resseguier | Geophysical flows under location uncertainty, Part I: random transport and general models[END_REF]. Hypothesis (2) can be removed, but the stochastic RTT expression assumes a more complex formulation. It is found that this constrain is naturally satisfied by fluids where density is constant in space (see section A.3), thus formula ( 24) is directly applied when the flow is incompressible or when the Boussinesq approximation is applied. On the contrary, an explicit formula has not been derived for a generic non-conserved quantity; hence hypothesis (3) is essential.

Notice that equation [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF] reduces to the classical RTT when the stochastic contribution in equation ( 1) is suppressed. This happens e.g. when σ ij = 0 and, consequently, the martingale η i t as well as the variance tensor are identically zero.

A concise derivation of the stochastic RTT is now presented. A generic random process φ(x, t) is expressed hereafter as a semimartingale of the form:

φ(x, t) = φ(x, 0) + t 0 g(x, s)ds + t 0 Ω f k (x, y, s)dB k s (y)dyds, (25)
where g, f are processes of bounded-variation and the Itō integral are used.

Proposition 1 (Differential of transported process) Let us consider φ a semimartingale of the type [START_REF] Métais | Spectral large-eddy simulation of isotropic and stably stratified turbulence[END_REF], sufficiently regular in space (bounded spatial gradient, two times derivable). If it is transported by a flow of the form [START_REF] Lesieur | Turbulence in Fluids[END_REF]; then, the time total-differential of φ is expressed by:

dφ(Xt, t) = ∂tφ + ∂φ ∂x i dX i t + 1 2 a ij (Xt, t) ∂ 2 φ ∂x i ∂x j dt + Ω σ ij (Xt, y, t) ∂ ∂x i f j (Xt, y, t)dydt, ( 26
)
where ∂t is the time partial-differential (i.e. with respect to second variable), and a ij is the variance tensor defined by equation ( 4).

Proof: The Itō -Wentzell formula is used to differentiate (in time) the transported process φ(Xt, t), corresponding to a composition of two processes. It reads:

dφ(Xt, t) = ∂tφ + ∂φ ∂x i dX i t + 1 2 d X i , X j t ∂ 2 φ ∂x i ∂x j + d ∂φ ∂x i , X i t , (27)
where the angular brackets denote the quadratic variation operation; e.g. see Le Gall [START_REF] Gall | Brownian Motion, Martingales, and Stochastic Calculus[END_REF] for an extended presentation. The following properties of the quadratic variation are recalled:

1. is symmetric and bilinear;

2. if g is a process of bounded-variation: g, B k t = 0

3. if f is deterministic function: f B i , B j t = f B i , B j t 4. singularity: d B i (y), B j (z) t = δ(y -z)δ ij dt
where Bt is a cylindrical Wiener process, δ(x) is the Dirac function and δ ij is the Kronecker symbol. Using these properties, the third and fourth terms in equation ( 27) are written explicitly. The third term is directly computed:

d X i , X j t = Ω σ ik (Xt, y, t)σ jk (Xt, y, t)dydt = a ij (Xt, t)dt. (28)
In the fourth term, the gradient of φ is obtained differentiating equation [START_REF] Métais | Spectral large-eddy simulation of isotropic and stably stratified turbulence[END_REF]:

∂ ∂x i φ(Xt, t) = ∂φ 0 ∂x i + t 0 ∂ ∂x i g(Xs, s)ds + t 0 Ω ∂ ∂x i f k (Xs, y, s) dB k s (y)dyds, (29)
then, the last term in ( 27) is rewritten as:

d ∂φ ∂x i , X i t = Ω σ ij (Xt, y, t) ∂ ∂x i f j (Xt, y, t)dydt. (30)
Substituting formula [START_REF] Moser | Direct numerical simulation of turbulent channel flow up to reτ =590[END_REF] and [START_REF] Neves | Wellposedness for stochastic continuity equations with ladyzhenskaya-prodi-serrin condition[END_REF] in equation ( 27), expression (26) is recovered.

Proposition 2 (Differential of transported and conserved process) Let us consider a stochastic process φ of the type [START_REF] Métais | Spectral large-eddy simulation of isotropic and stably stratified turbulence[END_REF]. If such a process is transported by a stochastic flow [START_REF] Lesieur | Turbulence in Fluids[END_REF] and is conserved, i.e. dφ(Xt, t) = 0, then:

∂tφ(Xt, t) = - ∂φ ∂x i w i dt + 1 2 a ij ∂ 2 φ ∂x i ∂x j dt + Ω ∂φ ∂x k ∂σ kj ∂x i σ ij dydt - Ω ∂φ ∂x i σ ik dB k t dy. ( 31
)
This formula expresses the time variation along a fluid-particle trajectory.

Proof: If φ(Xt, t) is conserved, then equation ( 26) can be re-arranged as follows:

∂tφ(Xt, t) = - ∂φ ∂x i w i dt - Ω ∂φ ∂x i σ ik (Xt, y, t)dB k t (y)dy- - 1 2 a ij ∂ 2 φ ∂x i ∂x j dt - Ω σ ij (Xt, y, t) ∂ ∂x i f j (Xt, y, t)dydt. ( 32
)
An expression of ∂tφ is obtained also from [START_REF] Métais | Spectral large-eddy simulation of isotropic and stably stratified turbulence[END_REF], and is compared with formula [START_REF] Oliveira | An improved piso algorithm for the computation of bouyancy driven flows[END_REF]. Exploiting the unique decomposition of the semimartingales, one obtains:

g(Xt, t) = - ∂φ ∂x i w i - 1 2 a ij ∂ 2 φ ∂x i ∂x j - Ω σ ij (Xt, y, t) ∂ ∂x i f j (Xt, y, t)dy, (33)
and the following implicit formula for f :

Ω f k (Xt, y, t) + ∂φ(Xt, t) ∂x i σ ik (Xt, y, t) dB k t (y)dy = 0. ( 34
)
This latter holds for every Brownian motion dB k t , thus:

f k (Xt, y, t) = - ∂φ(Xt, t) ∂x i σ ik (Xt, y, t). (35)
Substituting formula [START_REF] Pope | Turbulent Flows[END_REF] in equation ( 32) the final expression ( 31) is obtained.

Notice that a general form of a conserved semimartingale can be found by substituting equation ( 33) and [START_REF] Pope | Turbulent Flows[END_REF] in formula [START_REF] Métais | Spectral large-eddy simulation of isotropic and stably stratified turbulence[END_REF].

Translating equation [START_REF] Øksendal | Stochastic Differential Equations[END_REF] from Lagrangian to Eulerian coordinates and rearranging the second and third terms in the left-hand side, one obtains the expression of the material derivative (in differential form) within the stochastic framework.

Proposition 3 (Stochastic transport operator) If φ is a stochastic process of the type [START_REF] Métais | Spectral large-eddy simulation of isotropic and stably stratified turbulence[END_REF], transported by a stochastic flow [START_REF] Lesieur | Turbulence in Fluids[END_REF] and conserved, then the stochastic material derivative in differential form is:

Dφ(x, t) = ∂tφ + (w i dt + dη i t ) ∂φ ∂x i - 1 2 ∂ ∂x i a ij ∂φ ∂x j dt - 1 2
∂a ij ∂x i - Ω σ kj ∂σ ij ∂x i dy ∂φ ∂x j dt, (36)
which is reported in [START_REF] Resseguier | Geophysical flows under location uncertainty, Part I: random transport and general models[END_REF] as the stochastic transport operator for a conserved quantity.

The derivation of Stochastic RTT is now summarised. Let us consider a generic physical quantity, mathematically expressed by a stochastic scalar process q(x, t) that satisfies the hypotheses of Stochastic RTT. The solution of transport equation is found in the space of weak solutions.

Proof (Stochastic RTT): Consider a control volume V (t) and a test function ϕ(x, t) in the space domain Ω such that: it has compact support on V (t), it is conserved and satisfies [START_REF] Métais | Spectral large-eddy simulation of isotropic and stably stratified turbulence[END_REF]. Then, the weak transport equation for q reads: d

V (t) q(x, t)ϕ(x, t)dx = Ω ϕ∂tq + q∂tϕ + d q, ϕ t dx, (37)
applying the Itō integration by part and passing to the integral on Ω because ϕ has compact support on V . The last term on the right-hand side needs to be explicited. An expression of q and ϕ is given by the semimartingale decomposition [START_REF] Métais | Spectral large-eddy simulation of isotropic and stably stratified turbulence[END_REF]:

ϕ = ϕ(x, 0) + t 0 g(x, s)ds + t 0 Ω f j (x, y, s)dB j s dy, (38)
q = q(x, 0) + t 0 h(x, s)ds + t 0 Ω κ j (x, y, s)dB j s dy, (39)
where explicit formulae for g, h, f, κ are given, see proof of Proposition 2. Using these expressions to compute the quadratic variation, we get:

d q, ϕ t = d t 0 Ω κ i dB i dy, Ω f j dB j dz t = ∂q ∂x k ∂ϕ ∂x a k dt, (40)
The same expressions are differentiated to express ∂tq(x, t) and ∂tϕ(x, t), that are substituted in the transport equation ( 37) together with formula [START_REF] Resseguier | Stochastic modelling and diffusion modes for proper orthogonal decomposition models and small-scale flow analysis[END_REF]. Subsequently, ϕ is used to compute the weak derivative and gathered; the final equation reads:

d Ω q(x, t)ϕ(x, t)dx = Ω ϕ ∂tq + ∂qw i ∂x i dt + 1 2 ∂ 2 ∂x i ∂x j (qa ij ) dt- - ∂ ∂x k Ω q ∂σ kj ∂x i σ ij dydt + ∂ ∂x i qdη i t - ∂ ∂x ∂q ∂x k a k dt dx. (41)
Equation ( 41) is valid for every test function ϕ with compact support in V (t), thus:

d V q(x, t)dx = V ∂tq + ∂qw i ∂x i dt + ∂ ∂x i qdη i t + 1 2 q ∂σ ∂x 2 dt - 1 2 Ω q ∂σ ik ∂x j ∂σ kj ∂x i + ∂ 2 q ∂x i ∂x j σ ik σ kj + 2σ ij ∂q ∂x k ∂σ kj ∂x i dydt dx, (42)
where the terms are rearranged and the definition of variance tensor is used to simplify some terms. Equation ( 42) is the general form of Stochastic RTT, that is quite complex and eventually difficult to handle. One can notice that under the additional hypothesis that the random term is solenoidal in space, i.e.

∂ ∂x i dη i t (x) = 0 ⇔ ∂ ∂x i σ ik (x, y, t) ≡ 0, (43)
where the if and only if statement holds because the Brownian motion is arbitrarily chosen, the equation ( 42) simplifies to:

d V q(x, t) dx = V ∂tq + ∂ (qw i ) ∂x i dt - 1 2 ∂ 2 (qa ij ) ∂x i ∂x j dt + ∂ qdη i t ∂x i dx, (44)
that is the final form of Stochastic RTT.

In the following section, it is shown that the assumption of a solenoidal random turbulence field is satisfied by incompressible fluids, thus the simpler equation ( 44) can be used to derived the equation of motion.

It is worth noticing that the Stochastic RTT can be applied to all functions of the form [START_REF] Métais | Spectral large-eddy simulation of isotropic and stably stratified turbulence[END_REF]; specifically, to all process of bounded-variations that are a particular case of semimartingale where the martingale term is zero.

A.3 Derivation of Stochastic Navier-Stokes equations

The derivation of governing equations for fluid flows is performed with a similar strategy as in the classical framework, e.g. see [START_REF] Kundu | Fluid Mechanics[END_REF].

Conservation of mass

If ρ(x, t) is the mass density, then the conservation of mass is:

d V (t) ρdx = V (t) ∂tρ + ∂ρw i ∂x i dt + ∂ ∂x i ρdη i t + 1 2 ρ ∂σ ∂x 2 dt (45) - 1 2 Ω ρ ∂σ ik ∂x j ∂σ kj ∂x i + ∂ 2 ρ ∂x i ∂x j σ ik σ kj + 2σ ij ∂ρ ∂x k ∂σ kj ∂x i dydt dx = 0, (46)
where the general form of Stochastic RTT ( 42) is used here. For an incompressible fluid, density is constant ρ(x, t) = ρ and the mass conservation equation simplifies accordingly.

The integral is then removed exploiting the arbitrariness of control volume:

∂w i ∂x i dt + ∂ ∂x i dη i t + 1 2 ∂σ ∂x 2 dt - 1 2 Ω ∂σ ik ∂x j ∂σ kj ∂x i dydt = 0. ( 47
)
Separating the processes of bounded-variation and the martingales, the following system is recovered:

∂w i ∂x i + 1 2 ∂σ ∂x 2 - 1 2 Ω ∂σ ik ∂x j ∂σ kj ∂x i dy dt = 0, ∂ ∂x i dη i t = 0. (48)
Equation ( 48) shows that for an incompressible fluid the Brownian term is solenoidal; thus, the use of the simplified expression Stochastic RTT ( 44) is a posteriori justified for incompressible fluids. Using the solenoidal constraint, equations (48) leads to the system:

∂ ∂x i w i - 1 2 ∂ ∂x j a ij = 0, ∂ ∂x i σ ik = 0, (49)
which expresses the conservation of mass.

Conservation of momentum

Two derivations are proposed, they are named Lagrangian and Eulerian for convenience of notation. The former is based on the work of [27], the latter on that one of [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF].

• Lagrangian: The second Newton's law is:

d dt ρU i (Xt, t) = F i (Xt, t), (50)
where F i are the forces acting on a fluid-particle. If I i if the time integral of the forces (the impulse), equation ( 50) is re-written in a differential form as dρU i = dI i . It is expressed in a weak form as:

ρ h dU i (Xt, t) = h dI i (Xt, t), (51)
where h are test functions and ρ is constant. The left-hand side is:

ρ h dU i (Xt, t) = ρ h dw i (Xt, t) -ρ h dη i t , (52)
then, the right-hand side of equation (51) must have the same structure, see [27]. Hence, the impulse divides into two contributions and equation (51) becomes:

ρ h dw i (Xt, t) -ρ h dη i t = h dJ i -h dλ t t .

(53)

Matching similar terms, we arrived at the following relations:

ρdw i = dJ i , dη i t = dλ i t . (54)
Equations ( 54)-first is exploited to obtain the governing equation of motion, while (54)second states that the forces balance the contribution of random velocities. In Eulerial framework, this latter reads:

ρDw i (x, t) = dJ i (x, t), (55)
Applying the stochastic transport operator [START_REF] Pope | A lagrangian two-time probability density function equation for inhomogeneous turbulent flows[END_REF] with the solenoidal constrain (49), one gets:

Dw i = ∂tw i + (w j - 1 2
∂a kj ∂x k ) ∂w i ∂x j dt - 1 2 ∂ ∂x k a jk ∂w i ∂x j dt + dη j t ∂w i ∂x j , (56)
while the impulse is determined by a physical analysis of the forces acting on the system, as in the classical derivation:

dJ i = - ∂ ∂x i (pdt -
with µ is the fluid viscosity, and pressure is written in a semimartingale form [START_REF] Métais | Spectral large-eddy simulation of isotropic and stably stratified turbulence[END_REF] where dξt(x) = Ω ϑ i (x, y, t)dB i t (y)dy denotes the martingale contribution to pressure. Imposing the equality (55) and using the unique decomposition of semimartingale, the governing equations (63) are recovered.

• Eulerian: Once again, the momentum conservation is formulated in differential form. If

J i (x,
where the velocity decomposition (3) is employed. The impulse acting on V (t) is expressed by (58). Then, imposing equality (59) and separating the processes of bounded-variation to the martingales, one obtains: where the conservation of mass constrain (49) is applied to simplify the formula. Notice that the expected value of noise is zero because the random displacement is uncorrelated in time.

∂tw i + ∂w i w j ∂x j dt - 1
Then, the integral at the right-hand side can be interpreted as a spatial empirical mean of zero-mean random process, and have to be null. With this simplification, the system (63) of fluid dynamics equations is obtained.

Finally, the stochastic model for an incompressible (Newtonian) fluid reads: The system is composed by two coupled sets of deterministic and stochastic non-linear partial differential equations, in the unknowns w i and σ ij . The pseudo-stochastic model is obtained by avoiding the resolution of the last two stochastic equations, and closing the system by providing an expression a ij through physical assumptions.

                          
Let us also outline that the system (63) has been obtained under the assumption that the drift velocity is of bounded variation. Removing this assumption, the separation of the regular and the stochastic terms cannot be performed anymore. Hence, one obtains a fully stochastic Nevier-Stokes composed by stochastic partial differential equations. For geophysical flows (for isochoric flows in general), the continuity equations is also stochastic; see [START_REF] Resseguier | Geophysical flows under location uncertainty, Part I: random transport and general models[END_REF].

2 .

 2 The PSS does not rely on hypothesis (b). The tensorial form of a ij allows to reproduce the anisotropy of unresolved turbulence, i.e. different turbulent contributions along different directions. 3. The extra terms in PSS account for turbulent effects usually not considered in the classical models, namely turbulent advection and turbulent compressibility.

  ) turns into the classical solenoidal constrain. Therefore, the LES Smagorinsky model can be interpreted as a particular case of the PSS constant isotropic model.Approximation[START_REF] Bauer | Accuracy and robustness of nonlinear eddy viscosity models[END_REF] is valid if the turbulent energy is mainly concentrated in the region where the irrotational strain dominates vorticity. Exceptions on this behaviour have been found and have motivated the development of alternative models, like the Wall Adaptive Local Eddy-viscosity (WALE) model of Nicoud and Ducros[START_REF] Nicoud | Subgrid-scale stress modelling based on the square of the velocity gradient tensor[END_REF] or the structure function model of Métais and Lesieur[START_REF] Métais | Spectral large-eddy simulation of isotropic and stably stratified turbulence[END_REF].Approximation[START_REF] Brzeźniak | Stochastic partial differential equations and turbulence[END_REF] implies that the flow deformation rate can be represented by a linear function in each spatial point; thus it is a particularly regular function. This is equivalent to neglect the turbulent correction on advective velocity and continuity equation, hence the associated physical phenomena of turbophoresis and turbulent compressibility are not reproduced. 5 Numerical simulations PSS and LES are compared on turbulence channel flow at Re τ = 590. The former adopts a constant isotropic model for variance tensor, the latter adopts a constant Smagorinsky model for sub-grid scale viscosity. The Direct Numerical Simulation (DNS) of Moser et al. [28] is taken as reference.5.1 Case geometry and settingsThe channel is composed by two horizontal and parallel walls between which a shear flow develops. The dimensions in stream-wise (x), vertical (y) and span-wise (z) directions are 2πδ × δ × πδ, respectively. The flow is driven by a constant pressure gradient ∂p ∂x = -ρu τ /δ. The Reynolds number based on the friction velocity u τ is defined as Re τ = u τ δ/ν. The spatial variables are made non-dimensional as y + = yu τ /ν, the velocity as u + = u/u τ , time ast + = tu τ /ν.The characteristic flow time is estimated as t 0 = U 0 /2πδ, where U 0 is the bulk velocity in stream-wise direction.The computational domain is discretised by 96 × 96 × 96 points. They are uniformly distributed in stream-wise and span-wise directions, leading to a cell width ∆x + < 40 and ∆z + < 20, respectively. In vertical direction, the grid is stretched in a way such that the first cell is within y + = 1 and with 9 cells in y + ≤ 11; thus ensuring an accurate resolution of the boundary layer.

ing factor is set to λ = 5. 25 .

 25 Cyclic boundary conditions are set at the vertical boundaries, while velocity no-slip condition and pressure zero-gradient are imposed at the horizontal walls. All the cases are initialised with the instantaneous fields provided by a preliminary LES with constant Smagorinsky SGS model, that has reached the statistical steady state.5.2 Algorithm and implementationSimulations are performed taking advantage of the open-source software Open-FOAM v. 2.3.0. This is a C++ library for computational fluid dynamics and adopts the finite volume methods.The LESs are carried out using the solver pisoFoam included in the standard software distribution. The implementation details on this basic solver can be found in the official OpenFOAM documentation and in the work of Jasak et al.[15]. The constant Smagorinsky SGS model is provided by OpenFOAM, and its correct implementation was checked.Two PSSs are performed using the code pseudoStochasticPisoFoam, a home-made solver developed by the authors within the Fluminance research group at INRIA Rennes (France). The non-conservative form of pseudo-stochastic governing equations (5) are solved employing the Pressure-Implicit with Splitting of Operators (PISO) algorithm proposed by Issa et al.[START_REF] Issa | The computation of compressible and incompressible recirculating flows by a non-iterative implicit scheme[END_REF] and Oliveira

Figure 1 2 sFig. 1

 121 Figure 1 top-panel reports the mean non-dimensional stream-wise velocity along the wall coordinate. PSS and LES lead to similar profiles in the near-wall region (y + < 30), while the former exhibits slightly lower values in the log-law region (y + > 30). They underestimate the velocity magnitude at the centre channel and, as expected, both are not accurate in reproducing the boundary layer profile. This is a well known shortcoming of Smagorinsky model when c 2 s

Figure 1 Fig. 2

 12 Figure 1 bottom-panel displays the velocity RMS components. If ψ is a generic variable, we denote [ψ] rms = ψ 2 the root-mean square, where ψ = ψψ is the instantaneous fluctuation. Both PSS and LES collapse on the same profiles.Because the isotropic model is very similar to the Smagorinsky model, an improvement of accuracy by the PSS is not expected. The interest of this validation is to prove that the pseudo-stochastic model is as accurate as the state-of-the-art LES methodologies, despite its derivation relies to a substantially different framework and its governing equations include several extra terms, which are analysed in the following sections.

Figure 3 Fig. 3

 33 Figure 3 presents the non-dimensional turbulent advection velocity u + ta = u ta /u τ and the turbulent compressibility Φ + tc = Φ tc u 2 τ /ν are scrutinised. The stream-wise component of u ta is practically zero, as well as the span-wise component; thus they are not displayed. The vertical component profile reveals low

Figure 4

 4 Figure4displays the Φ tc negative (blue) and positive (orange) isosurfaces near the bottom wall, at an instantaneous flow configuration. They are organised in spots, confined in the near-wall region and elongated in the streamwise direction. In accordance with the Φ tc mean profile, the negative spots are closer to the wall (y + < 10), while the positive one are immediately above (10 < y + < 20). The shape and the location of the isosurfaces suggest a correlation with the streaks structures that characterises turbulent wall flows. The streaks are generated in a region of low velocity, very close to the wall, approximately at y + 5. They are elongated in the stream-wise direction, with a characteristic length of ∆x + ∼ = 1000 and a span-wise period of ∆z + ∼ = 100. This estimation can vary with respect to the wall distance, see Smith & Metzler[START_REF] Smith | The characteristics of low-speed streaks in the near-wall region of a turbulent boundary layer[END_REF]. Despite their widespread presence, there is no clear consensus on the streak formation mechanism and multiple theories have been proposed in literature, see Chernyshenko & Baig[START_REF] Chernyshenko | The mechanism of streak formation in near-wall turbulence[END_REF]. The Φ tc isosurfaces have, overall, the

Fig. 4 Fig. 5

 45 Fig. 4 Positive and negative isosurfaces of Φtc near the bottom wall at an instantaneous flow configuration. Orange: isosurface at Φ + tc = 3.5/times10 -4. Blue: isosurface at Φ + tc = -3.5/times10 -4.

  the classical large-eddy simulations (LES) approach. The PSS model is based on an innovative decomposition of the fluid-particle trajectory in a drift displacement and a stochastic perturbation. The former reproduces the mean flow, the later accounts for the turbulent perturbations which are modeled as a Brownian motion. Imposing such a decomposition, together with a regularity assumption on the drift velocity, a set of deterministic and stochastic equations of motion are derived using stochastic calculus; then, the pseudostochastic equations are obtained by neglecting the solution of stochastic equations and closing the system by physical assumptions. The result is a new set of governing equations which includes extra terms deriving from the stochastic modeling of turbulence. The PSS model is found to be a generalisation of the classical Navier-Stokes equations, and reproduces phenomena usually not considered: turbophoresis and turbulent compressibility.

  t) is the impulse of total forces per volume, the second law of mechanics reads:

			d	ρU i (x, t)dx =		dJ i (x, t)dx.	(59)
				V (t)				V (t)
	The Stochastic RTT is applied to the left-hand side in order to get:
	d	V (t)	ρU i dx =	V (t)	ρ ∂tw i +	∂ ∂x j	(w i w j )dt -	1 2	∂ 2 ∂xs∂x k	(w i a sk ) dt+
			+	∂ ∂x j	w i dη j t	dx + d	V (t)	ρ ηi t dx,