Carlo Cintolesi
email: carlo.cintolesi@gmail.com

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

Pseudo-stochastic simulation of turbulent channel flows with near-wall modelling

Keywords: Stochastic model, Turbulence, Near-wall models, Numerical simulations, OpenFOAM

) is investigated and compared with the large-eddy simulation methodology. The theoretical analysis shows that this model is a generalisation of the eddy-viscosity model, which does not undergo the same restrictive physical assumptions and describes physical phenomena usually not considered (turbophoresis and turbulent compressibility). Numerical simulations of turbulent channel flows are performed. In order to better reproduce the turbulence anisotropy, a near-wall damping function is derived and successfully validated: the damping is imposed only on wall-normal direction (minimal constraint) and it requires to set a single parameter (reduced empirical content). Simulations show the accuracy of the new model, especially when the computational grid becomes coarse. A weak turbophoresis phenomenon is detected near the wall, while turbulent compressibility effects appear to be possibly related to the streaks structures.

Introduction

The use of stochastic calculus to describe fluid flows appears to be a suitable strategy for turbulence modelling in computational fluid dynamics. The random nature of turbulence cannot be completely represented by means of deterministic variables, while it is the specific purpose of stochastic processes. Nevertheless, the numerical solution of stochastic equations and the mathematical complexity inherent to the use of stochastic calculus poses challenging issues. Turbulence modelling with stochastic variables is of great interest in geophysical flow analysis, where the unresolved processes related to coarse spatial discretisation are handled with probabilistic models. In the same spirit, stochastic models can be applied to numerical simulations of environmental and engineering flows.

In the last decades, several efforts have been made in this concern. In the context of the Probability Density Function (PDF), the Langevin equation was used to describe the velocity of a fluid particle subject to a turbulent flow, modelled as a Brownian motion (see [START_REF] Pope | Turbulent Flows[END_REF]). First applications focused on homogeneous isotropic turbulence; later extended by [START_REF] Pope | A lagrangian two-time probability density function equation for inhomogeneous turbulent flows[END_REF] to inhomogeneous case and by [START_REF] Durbin | Realizability of second-moment closure via stochastic analysis[END_REF] to anisotropic diffusion case. In the Eddy-Damped Quasi-Normal Markovian (EDQNM) models, introduced by [START_REF] Orszag | Analytical theories of turbulence[END_REF] and [START_REF] Leslie | Developments in the theory of turbulence[END_REF], the large-scale governing equations were closed in the spectral space by modelling the third/fourth-order moments through a Gaussian closure. This strategy was found to be suitable in case of strong non-linearity in the small-scale turbulence. [START_REF] Chasnov | Simulation of the kolmogorov inertial subrange using an improved subgrid model[END_REF] developed a forced-dissipative model, where the large-eddy Navier-Stokes equations were corrected by an eddy-viscosity and a stochastic force terms. Similarly, [START_REF] Leith | Stochastic backscatter in a subgridscale model: Plane shear mixing layer[END_REF] studied the case of plane shear mixing layer and improved the accuracy of LES with the Smagorinsky model by adding an empirical stochastic backscatter. The work of [START_REF] Kraichnan | Dynamics of nonlinear stochastic systems[END_REF] exploited a different approach: the Navier-Stokes equations were replaced by a set of equations having the same mathematical properties, closed by a Gaussian stochastic model. This model led to valuable results when applied to the study of mathematical properties and physical effects, like turbulent diffusion and backscatter. [START_REF] Frederiksen | Subgrid modelling for geophysical flows[END_REF] showed that the same methodology can be used in the stochastic modelling of barotropic flows or in quasi-geostrophic approximation, as well as for the description of the interactions between topography and small-scale eddies.

Such attempts to include random functions in fluid dynamics modelling exhibit some limitations: in POF and EDQNM models the solution is found in the spectral space instead of the physical one; the explicit introduction of random term relies mostly on empirical considerations and leads to a certain degree of arbitrariness. For example, a question arises whether the random forcing term should be multiplicative or additive.

An alternative approach was developed. It is based on the idea that the velocity field itself is a random process, composed of a differentiable component and a fast oscillating random term. Physically, the former describes the smooth macroscopic velocity while the latter accounts for the stochastic turbulent motion. Under this assumption, the fluid dynamics equations are re-derived using stochastic calculus, leading to a complete set of stochastic partial differential equations. Pioneering work in this sense was made by [START_REF] Brzeźniak | Stochastic partial differential equations and turbulence[END_REF].

Subsequently, [START_REF] Mikulevicius | Stochastic navier-stokes equations for turbulent flows[END_REF] and [START_REF] Flandoli | The interaction between noise and transport mechanisms in pdes[END_REF] expanded his formulation and studied the mathematical properties of the resulting stochastic system. Such a model has been further developed by [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF] in view of practical applications and takes the name of model under Location Uncertainty (LU). Later, [START_REF] Neves | Wellposedness for stochastic continuity equations with ladyzhenskaya-prodi-serrin condition[END_REF] theoretically investigate a similar system, while [START_REF] Holm | Variational principles for stochastic fluid dynamics[END_REF] derives an equivalent model using Lagrangian mechanics. This last model differs from LU because an extra term appears in the momentum equation, which ensures helicity and circulation conservation but may alter the kinetic energy budget.

The LU model was tested in several cases: Resseguier et al. (2017a,b,c) successfully used such type of model to study geophysical flows, which was found to be more accurate in the reproduction of extreme events and provided new analysis tools. [START_REF] Chapron | Largescale flows under location uncertainty: a consistent stochastic framework[END_REF] investigated the Lorentz-63 case and state that LU is more effective in exploring the regions of the deterministic attractor than the classical models. Furthermore, it was used in conjunction with the proper orthogonal decomposition technique by Resseguier et al. (2017d) for studying a wake flow past a circular cylinder at Re = 3900. Recently, [START_REF] Pinier | A model under location uncertainty for the mean velocity in wall bounded flows[END_REF] perform mathematical analysis of the turbulent boundary layer through the LU equations. They propose a complete explicit profile for the mean vertical velocity, that includes an expression for the velocity in the buffer layer, for which a rigorous theoretical model is missing so far.

Despite these encouraging results, to perform stochastic numerical simulations for practical applications poses some difficulties; e.g. the numerical resolution techniques are not straightforward and they can possibly require a large computational effort. In order to circumvent such difficulties, [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF] introduced a hybrid model hereafter named pseudo-stochastic model : first, the governing equations are decomposed into two coupled system of partial and stochastic differential equations; second, the resolution of the latter is avoided and the system is closed by modelling the effects of the random velocity term through physical assumptions. Hence, the flow dynamics is described by a set of classical partial differential equations, which includes terms that derive from the stochastic representation of turbulence. [START_REF] Harouna | Stochastic representation of the reynolds transport theorem: Revisiting large-scale modeling[END_REF] used the pseudo-stochastic model to investigate the Green-Taylor vortex flow, testing several closure models. [START_REF] Chandramouli | Coarse large-eddy simulations in a transitional wake flow with flow models under location uncertainty[END_REF] successfully employed the model to simulate the transitional wake flow with coarse mesh resolution.

The present contribution aims to explore the potentiality of the pseudostochastic model, making a direct comparison with the Large-Eddy Simulation (LES) methodology. First, the model is described and discussed in details; then numerical simulations on the turbulent channel flows are performed and analysed. The main novelties here reported are: a detailed study of the pseudostochastic equations with respect to the classical ones; the derivation of a (resolved) turbulent kinetic energy budget for LU; the development of a near-wall model for pseudo-stochastic simulations from the study carried out by [START_REF] Pinier | A model under location uncertainty for the mean velocity in wall bounded flows[END_REF].

The paper is organized as follows: section 2 describes the pseudo-stochastic equations, along with the kinetic energy budget and the near-wall model; section 3 reports a physical interpretation of the equation terms, as well as a comparison with the LES methodology; section 4 presents the simulation methodologies and settings; section 5 discusses the validation of the near-wall model and the simulation results; section 6 reports some final remarks.

Pseudo-stochastic model

The pseudo-stochastic equations are described, together with the kinetic energy budgets. We refer to [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF] and [START_REF] Resseguier | Mixing and fluid dynamics under location uncertainty[END_REF] for the formal derivation.

Stochastic formalism

The pathlines in a turbulent flow are modelled as a stochastic process, where a regular function is perturbed by a random (turbulent) process. Consequently, a Lagrangian fluid-particle displacement is described by a stochastic differential equation of the type:

dX i t (x 0 ) = w i (X t , t)dt + Ω σ ik (X t , y, t)dB k t (y) dy, (1)
where the index i = 1, 2, 3 indicates respectively the x,y,z-component in the space domain Ω (they are placed at top or bottom indifferently) and the Einstein summation convention is adopted; X i t is the trajectory followed by a fluidparticle initially located in x 0 ; w i is a differentiable function that corresponds to the drift velocity; dη i t = Ω σ ik dB k t dy is a stochastic process (accounting for turbulent effects) uncorrelated in time but correlated in space. This last is constructed as a combination of a cylindrical Wiener processes B k t (x) not differentiable in time, and a time-differentiable symmetric diffusion tensor σ ik (x, y, t)

which acts as an integral kernel. Hence, they are fast oscillating stochastic components, possibly anisotropic and inhomogeneous in space.

The velocity field U i in Eulerian coordinate x is derived from equation (1):

U i (x, t) = w i (x, t) + ηi t (x), (2)
where the second term on the right-hand side expresses the stochastic velocity defined as the weak derivative of η i t (x) in time. 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 (2003).

In the derivation of the stochastic model, the quadratic variation of the diffusion tensor is of particular interest since it represents the time-variation of spatial variance of the stochastic increments along time. It is named as the variance tensor and it is defined as:

a ij (x, t) = Ω σ ik (x, y, t)σ kj (x, y, t) dy, (3)
it can be shown to be a symmetric and semi-positive definite matrix with dimension [m 2 /s].

Pseudo-stochastic equations of motion

The stochastic process (1) that described the flow is not time-differentiable in the framework of classical analysis. Thus, the Navier-Stokes equations need to be re-derived using the stochastic calculus, where the use of the Itō-Wentzell formula is crucial for computing the derivative in time, see [START_REF] Kunita | Stochastic Flows and Stochastic Differential Equations[END_REF]. The result is a complete system of stochastic partial differential equations that describes the fluid flow. Assuming the drift velocity is of bounded variation (deterministic) and using the unique decomposition of semi-martingale, the system can be divided into a set of stochastic equations and a set of pure deterministic ones. The former allows finding an expression for the variance tensor a ij , required for the resolution of the latter. The pseudo-stochastic model is derived by neglecting the resolution of the stochastic equations and closing the system by giving an expression of the variance tensor, which is modelled through physical hypothesis. This choice gives rise to a hybrid model where the terms that depend on a ij accounts for the Stochastic Unresolved Scales (SUS) of motion.

The pseudo-stochastic equations for incompressible flows read:

       ∂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.
(4) they represent the momentum and mass conservation, respectively, written in the non-conservative form proposed by Resseguier et al. (2017a). The effective advection velocity w * is defined as:

w * i = w i - 1 2 ∂a ik ∂x k , (5)
and the pressure is the sum of an hydrostatic pressure and an isotropic turbulent term:

p = p h + ν 3 ∂w ∂x = p h + ν 6 ∂ 2 a sk ∂x k ∂x s . ( 6
)
This last 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].

It is worthwhile to notice that system (4) reduces to the classical Navier-Stokes equations when the variance tensor tends to the zero matrix, i.e. when the stochastic contributions disappear.

In the framework of computational fluid dynamics, the drift velocity w i can be interpreted as the (numerically) resolved velocity field, while the random field η i t assembles the (turbulent) unresolved motions. Therefore, giving an expression on variance tensor is equivalent to specifying a turbulence model.

Resolved kinetic energy budget

Equations for mean and turbulent kinetic energy budget of the resolved scales of motion are here derived. The resolved velocity is decomposed in a mean and a fluctuating part, respectively:

w i = W i + w i , (7)
where the capital 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 variance tensor accounts for the SUS effects on the mean flow.

The budget of resolved kinetic energy K = (W i W i )/2 is obtained multiplying momentum equation (4-first) by W i and averaging. Applying the conservation of mass (4-second) and rearranging the terms, one gets:

∂K ∂t + W j - ∂ ∂x k A jk 2 ∂K ∂x j = ∂ ∂x j -P W j + (νδ jk + A jk 2 ) ∂K ∂x k (8) -w j - ∂ ∂x k a jk 2 w i W i + a jk ∂w i ∂x k W i (9) + p 2 ∂ 2 A jk ∂x j ∂x k -νδ jk + A jk 2 ∂W i ∂x j ∂W i ∂x k (10) + w j - ∂ ∂x k a jk 2 w i ∂W i ∂x j - a jk 2 ∂w i ∂x j ∂W i ∂x k (11)
The second term on the left-hand side represents the rate of change by means of the effective (mean) advection. The first four terms on the right-hand side express the energy transport by pressure, molecular and turbulent viscous stresses, resolved turbulence, turbulent SUS motion (respectively). The fifth term is due to the non-solenoidal velocity field and is related to the compression-expansion work made by the SUS; it can be a production or dissipation term. The sixth term is a viscous and turbulent dissipation (it can be proven that A ij is positive defined), while the seventh term is a loss due to resolved turbulence; the same term but with opposite sign is present in the turbulent kinetic energy budget presented later in this section. The last term indicates dissipation/production due to SUS.

The (resolved) turbulent kinetic energy κ = w i w i /2 budget is obtained following the procedure described in [START_REF] Kundu | Fluid Mechanics[END_REF]: the equation for resolved fluctuations is obtained subtracting expression ( 8) from (4-first), then multiplying by w i and averaging. Using the continuity equation (4-second) to simplify the terms and rearranging them, one obtains the following expression for stochastic Turbulent Kinetic Energy (TKE):

∂ κ ∂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 (12)
On the left-hand side, the second and third terms represent the TKE advection by mean and SUS effective advection velocity, respectively. On the right-hand side:

• the first four terms express spatial transport;

• the fifth term is a turbulent compression/expansion term due to SUS;

• the sixth and seventh terms account for dissipation by molecular viscosity, resolved turbulence and SUS motions;

• the eight term represents the shear production, including the contribution by the fluctuations of turbulent advection velocity;

• the last term indicates a loss due to SUS; this term is also present in the resolved kinetic energy budget.

Both the kinetic energy and TKE expressions reduce to the classical ones if the stochastic contribution is negligible a ij 0.

Isotropic constant model for variance tensor

Several strategies can be adopted to model the variance tensor. The isotropic model is developed by analogy with the Smagorinsky model, e.g. see [START_REF] Deardorff | A numerical study of three-dimensional turbulent channel flow at large reynolds numbers[END_REF], and was first proposed by [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF]. The variance tensor is given by:

a ij = c m ∆ 2 |S| δ ij , (13)
where c m is a model parameter, |S| is the strain-rate tensor norm, and ∆ is the computational cell width. The variance tensor reduces to a diagonal matrix with equal elements because turbulence is assumed isotropic and homogeneous in all directions.

Near-wall modelling of variance tensor

In a very recent work, [START_REF] Pinier | A model under location uncertainty for the mean velocity in wall bounded flows[END_REF] studied the mean velocity profile of the turbulent boundary layer through the LU equations. They proposed a modification of the classical velocity expression for wall-bounded flow and provided an analytical formula for the buffer layer, not available till now. Notice that the modified advection velocity plays a crucial role in the mathematical derivation of this formula; therefore, such a profile cannot be deduced using the classical formulation of the Navier-Stokes equations, where the modified advection is not explicitly taken into account. In the viscous sublayer (y

+ < y + 0 )
and in the logarithm region (y + L < y + < y + 1 ) the linear and log-law velocity profiles (respectively) are retrieved, while in the buffer layer (y + 0 < y + < y + L ) a hyperbolic profile is specified:

u + (y + ) =                y + y + ∈ [0, y + 0 ] u + (y + 0 ) + 2 κ - 4 κ ( κy + -y + 0 ) + 2 y + ∈ [y + 0 , y + L ] u + (y + L ) + 4y + L κ(y + L -y + 0 ) + 2 2 ln y + y + L y + ∈ [y + L , y + 1 ] (14)
where u(y) is the streamwise velocity as a function of the wall-normal coordinate.

Quantities are made non-dimensional by means of the friction velocity u τ and molecular viscosity ν, as usual: y + = yu τ /ν and u + = u/u τ . The κ is a model constant (to not be confused with the von Kármán constant); for a plain channel flow it has been estimated to be κ = 0.158 from direct numerical simulations.

The boundaries of the three regions are: An additional result concerns the expression of the variance tensor. In the viscous sublayer, a ij is almost zero, while in the buffer layer the wall-normal component depends only from the distance from the wall and exhibits a linear profile:

a + yy (y + ) = κ y + -y + 0 , (15)
where a + ij = a ij /ν. In the log-law region, it scales as the square-root of the wall distance:

a + yy (y + ) = κ y + L -y + 0 y + /y + L . (16)
No estimations are provided for the other components.

Preliminary pseudo-stochastic simulations with the isotropic constant model ( 13) have shown an excessive energy dissipation near the solid boundaries, given by high values of a ij in the buffer and viscous layer. This is not unexpected since the LES Smagorinsky model (that is the classical counterpart of the isotropic model) exhibit the same behaviour (see discussion in following section 3.2).

To correct this behaviour, a damping function for variance tensor is here formulated, exploiting the above-described characterisation of wall-normal component. Away from the wall, a yy is given by the isotropic model; at a point y + B placed in the buffer layer, a linear decrease is imposed in such a way to reach the zero value at y + 0 ; in the viscous sublayer, it is set to be zero. Hence, the LU near-wall model reads:

a + yy (y + ) =              0 y + ∈ [0, y + 0 ] a + yy (y B ) y + -y + 0 y + B -y + 0 y + ∈ [y + 0 , y + B ] c m ∆ 2 ν |S| δ ij y + ∈ [y + B , y + 1 ] (17)
The coordinate y + B is a model parameter that have to be set after theoretical or numerical estimation. No constraints are imposed on the other components; they are computed according to the isotropic model (13).

Physical interpretation and comparison with LES models

The pseudo-stochastic equations (4) are analysed from a physical point of view, and a comparison with the eddy-viscosity model used in LES is reported.

Physical interpretation

Recalling the decomposition of the velocity gradient as the sum of the symmetric and the antisymmetric part, respectively the strain-rate tensor S ij = 1 2 (∂w i /∂x j + ∂w j /∂x i ) and the rotation-rate tensor Ω ij = 1 2 (∂w i /∂x j -∂w j /∂x i ), the pseudo-stochastic equations (4) are rearranged as:

∂w i ∂t + w j - 1 2 ∂a jk ∂x k ∂w i ∂x j = - ∂p ∂x i + ∂ ∂x j 2νδ jk + a jk 2 S ki - ∂ ∂x j a jk 2 Ω ki , (18)
and

∂w i ∂x i = 1 2 ∂ 2 a jk ∂x j ∂x k . ( 19
)
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. As pointed out by Resseguier et al. (2017a), it corresponds to a velocity induced by the unresolved eddies, that is 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.

Diffusion due to SUS: the last two terms on the right-hand side of equation ( 18) account for the turbulent diffusion; the variance tensor plays the role of a diffusion tensor similar to a generalised eddy-viscosity coefficient.

Both the deformation rate and rotation-rate contribute to diffusion, unlike in the classical eddy-viscosity model in which fluid rotation-rate is assumed to be irrelevant in turbulent modelling.

Turbulent compressibility: the continuity equation ( 19) suggests that the flow is turbulent-compressible; i.e. the unresolved turbulence induces a local fluid compression or expansion.

The variance tensor (3) 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 of the SUS. The role played in governing equations ( 4)

and in kinetic energy budgets ( 8)-( 12), suggests that 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 (increasing or possibly decreasing) of fluid viscosity, which 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 ; ( 20
)
the divergence of the turbulent advection velocity measures the turbulent compressibility:

Φ tc = 1 2 ∂ 2 a ij ∂x i ∂x j , ( 21
)
and it is directly proportional to the isotropic turbulent term appearing in equation (6).

Comparison with LES eddy-viscosity models

In the classical framework, the fluid velocity u(x, t) is a deterministic function of time and space. Adopting the LES approach, the computational grid act on the governing equations as an implicit spatial filter (denoted by an over-bar)

depending on the local cell width ∆ = (∆x∆y∆z) 1/3 , see [START_REF] Sagaut | Large eddy simulation for incompressible flows[END_REF] and [START_REF] Piomelli | Large-eddy and direct simulation of turbulent flows[END_REF] for extended reviews. Filtering the Navier-Stokes equations, the sub-grid scale (SGS) stress tensor τ ij = (u i u ju i u j ) appears:

       ∂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 (22)
and it has to be modelled to close the system: a popular choice is to use the eddy-viscosity models. They are a class of turbulent models relying on the Boussinesq assumption, where the anisotropic part of τ ij is proportional to the resolved strain-rate tensor through ν sgs the SGS viscosity parameter:

τ R ij = τ ij - τ kk 3 δ ij = -2ν sgs S ij , (23)
while the isotropic part is included in the pressure term and does not contribute to the motion. This parameter has to be specified by additional models; the classical constant Smagorisky model is here analysed:

ν sgs = c 2 s ∆ 2 S , (24)
where S is the norm of the filtered strain-rate tensor, and the parameter c 2 s is set constant and can be evaluated from experiments, direct numerical simulations or analytical considerations, e.g. see [START_REF] Lilly | The representation of small-scale tubulence in numerical simulation experiments[END_REF]. The main drawback of this approach is to rely on the homogeneous turbulence assumption. This hypothesis is violated in many, even simple, cases. For example, close to solid surfaces where the turbulent length-scales decrease. To cope with this shortcoming, a damping function is usually introduced in order to account for the reduction of turbulence intensity. After the first work of van Driest (1956) several modifications of the original damping function have been proposed, e.g. see Piomelli et al. (1989) Cabot and[START_REF] Cabot | Approximate wall boundary conditions in the largeeddy simulation of high reynolds number flow[END_REF]. They can be summarised in the following expression:

∆ = min κy C δ 1 -e -y + A + n m , ∆ , (25)
where κ = 0.41 is the von Kármán constant. The original formulation by van [START_REF] Van Driest | On turbulent flow near a wall[END_REF] prescribe n = m = 1, A + = 0.26 and C δ = 1.00.

Remarks on eddy-viscosity model

Notice that the eddy-viscosity equation ( 23) implies that the Boussinesq's hypotheses are satisfied: (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 an eddy-viscosity model if the variance tensor is expressed by a ij = 2ν sus δ ij , i.e. assuming that the SUS induce an (isotropic) increasing of fluid viscosity. In this sense, the pseudo-stochastic model can be considered as a generalisation of the eddy-viscosity model. The comparison between the two models points out some theoretical advantages of the former:

(i) The effects of unresolved scales of motion are given by a ij , without imposing any constraints on the directions along with the SUS acts on the resolved flow. Hence, hypothesis (a) is not required.

(ii) The tensor form of a ij allows reproducing the anisotropy of unresolved turbulence, i.e. different turbulent contributions along different directions.

Thus, hypothesis (b) is not required.

(iii) The extra terms in the governing equations account for turbulent effects usually not considered in the classical models, namely turbulent advection and turbulent compressibility.

The eddy-viscosity models are 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 [START_REF] Pope | Turbulent Flows[END_REF] for an overview on this issue.

It is worth mentioning that the eddy-viscosity parameter a ij comes directly from the basic assumption of the velocity decomposition (2); whereas it is introduced in LES equations through an ad hoc physical assumption. Overall, the pseudo-stochastic model represents a general approach that overcomes the limitations of the Boussinesq assumption and includes turbulent effects not considered in the classical LES sub-grid scales models.

Remarks on Smagorinsky model

Expanding the pseudo-stochastic isotropic model ( 13), it can be shown that it reduces to the LES Smagorinsky model under two approximations:

(i) the rotation-rate does not contribute to turbulence effects on the mean flow;

(ii) the norm of the strain-rate tensor is almost harmonic (Laplacian is close to zero).

Notice that with the latter hypothesis the continuity equation (4-second) boils down to the classical solenoidal constraint. Therefore, the LES Smagorinsky model can be interpreted as a particular case of the pseudo-stochastic isotropic model.

Approximation (i) 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 [START_REF] Nicoud | Subgrid-scale stress modelling based on the square of the velocity gradient tensor[END_REF] or the structure function model of [START_REF] Métais | Spectral large-eddy simulation of isotropic and stably stratified turbulence[END_REF].

Approximation (ii) 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.

Simulation methodologies

The LU near-wall model ( 17 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, numerical solvers are secondorder accurate in time and space. The time advancement fulfils the Courant-Friedrichs-Lewy condition Co < 0.5. For LES, the Courant number is computed as Co = ∆t|u|/δx, where ∆t is the time step, |u| is the velocity magnitude through the cell, δx is the cell length. For PSS, the definition is modified in order to account for the effective advection velocity: Co = ∆t|w * |/δx.

Case geometry and settings

The channel is composed of two horizontal and parallel walls between which a shear flow develops. The dimensions in stream-wise (x), vertical (y) and spanwise (z) directions are 2πδ × δ × πδ, respectively. Several discretisation meshes are employed, whose parameters are summarised in Table 1. The computational points are uniformly distributed in streamwise and spanwise directions, while the grid is stretched along the vertical direction. The stretching is symmetric with respect to the channel center plane y = δ, and it is obtained with a doubleside stretching function based on hyperbolic tangent:

y(ξ) = 1 2 1 + tanh(λ(ξ -1/2)) tanh(λ/2) , ( 26
)
where ξ is the vertical coordinate of uniform point distribution. [START_REF] Moser | Direct numerical simulation of turbulent channel flow up to re τ =590[END_REF]. Top profiles: dash violet, analytical profile (17) derived from LU by [START_REF] Pinier | A model under location uncertainty for the mean velocity in wall bounded flows[END_REF]. Bottom profiles: red symbols, PSS with near-wall model; red lines, LES with van Driest damping; blue symbols, PSS without near-wall model; blue lines, LES without van Driest damping.

< u + x > y + ÅÍÄÍ¸Ú ×ÓÙ× ÅÍÄ͸ Ù Ö ÅÍÄ͸ÐÓ ¹Ð Û AEË ÈËË Û ÐÐ ÑÓ Ð Ä Ë Ú Ò Ö ×Ø ÈËË ÒÓ ÑÔ Ò Ä Ë ÒÓ ÑÔ Ò

Non-dimensional parameters

Quantities are made non-dimensional by the friction velocity u τ and molecular viscosity ν as follow: space x + = xu τ /ν, time t + = tu 2 τ /ν, velocity

u + = u/u τ , variance tensor a + ij = a ij /ν.
The flow is driven by a constant pressure gradient ∂p ∂x = -ρu τ /δ; Reynolds number is set to Re τ = u τ δ/ν. The characteristic flow time is estimated to be

t 0 = U 0 /2πδ
, where U 0 is the bulk velocity in stream-wise direction, while the large-eddy turn over time is estimated to be t * = tu τ /δ.

Results and discussion

The following notation is adopted: if φ is a generic variable, then φ is the time and space averaged over x, z-directions, φ = φφ is the instantaneous fluctuation and [φ] rms = φ 2 is the root-mean square. After the statistical steady state is reached, statistics are collected in an interval of 30t * ∼ 3t 0 every 0.1t * .

Near-wall model assessment

The LU near-wall model for variance tensor is validated in the plane channel flow Re τ = 395. The computational grid is described in Table 1, and ensures a high resolution of the flow. Four simulations are performed: PSS that enforce the near-wall model, LES with van Driest damping, PSS and LES switching off the near-wall models. Figure 1 shows the non-dimensional mean streamwise velocity along vertical direction. In the top-plot, analytical expression ( 14) for mean velocity is compared with the DNS data: in all the three boundary layer regions, the velocity profile is correctly described. Particularly, there is a good agreement between the hyperbolic function and the reference data in the buffer layer. In the bottom-plot, the results of the simulations with and without near-wall models are reported. The data of the PSS and LES collapse one onto the other; hence they are discussed together. As expected, when the near-wall models are disabled, the velocity profile is underestimated. This is caused by a non-physical high level of eddy-viscosity near the wall (see also discussion of Figure 2), that induces a large energy dissipation. When the near-wall models are activated, velocity is well captured in the viscous and buffer layer.

Figure 2 presents the non-dimensional mean eddy-viscosity parameters for LES and PSS, respectively ν + sgs = ν sgs /ν and a + ij = a ij /ν. Simulations are compared with the SGS eddy-viscosity profile reported in [START_REF] Armenio | A lagrangian mixed subgrid-scale model in generalized coordinates[END_REF]. Such a profile is obtained from LES of the channel at Re τ = 395 with a spectral code described in [START_REF] Sarghini | Scale-similar models for large-eddy simulations[END_REF]. The size and the discretisation of the computational domain are comparable to the one used here. The spectral code implements the Lagrangian dynamic model of [START_REF] Meneveau | A lagrangian dynamic subgrid-scale model of turbulence[END_REF],

where the eddy-viscosity is computed cell-by-cell, by comparing two scales of

< u + x > y + ÈËË Ú ÖÝ Ó Ö× Ä Ë Ú ÖÝ Ó Ö× ÈËË Ó Ö× Ä Ë Ó Ö× ÈËË Ò Ä Ë Ò AEË
Figure 3: Non-dimensional mean velocity profiles along wall-normal direction of turbulent channel at Reτ = 590, for the three meshes described in Table 1. Solid lines, DNS by [START_REF] Moser | Direct numerical simulation of turbulent channel flow up to re τ =590[END_REF]; symbols, PSS; dash lines, LES.

Overall, the PSS with the near-wall model is able to reproduce the velocity profile as well as the LES. The eddy-viscosity profile is correctly reproduced near the wall, where the LU damping is in a good agreement with reference data. Notice that the damping is imposed only on the wall-normal component of the generalised eddy-viscosity tensor (a yy ), while no modification are required for the wall-parallel components a xx , a zz . Moreover, the only parameter to be set is the damping point y + B . On the contrary, the class of van Driest functions (25) are applied to all the velocity components and required to choose several empirical parameters, which leads to larger empirical content.

Channel flow analysis

The PSS with near-wall damping is compared with LES van Driest damping on three different meshes with a decreasing resolution in wall-parallel directions (see Table 1) at Re τ = 590.

Figure 3 displays the non-dimensional streamwise velocity component. PSS and LES practically collapse on the same values. They exhibit accurate results in the inner-region (y + < 50) for all the meshes; whereas they tend to overestimate velocity in the outer-region (50 < y + ). Such overestimation increases as the computational grid degrades. For a very coarse grid, the PSS shows a slightly better profile with respect to LES in the buffer layer (10 < y + < 30), as a consequent of a different damping (see Figure 5)

Figure 4 reports the root-mean square (RMS) of velocity components. In general, no significant differences are detectable between PSS and LES. As expected, the profiles are more accurate as the mesh resolution increases. The streamwise RMS is overestimated and the peak moves from the buffer layer towards the log-law region as the mesh becomes coarser and coarser. Notice that

[u + z ] rms y + Ä Ë Ò ÈËË Ò AEË [u + y ] rms Ä Ë Ó Ö× ÈËË Ó Ö× [u + x ] rms Ä Ë Ú ÖÝ Ó Ö× ÈËË Ú ÖÝ Ó Ö×
Figure 4: Non-dimensional root-mean square of velocity components along wall-normal direction of turbulent channel at Reτ = 590. Simulations with the three meshes described in Table 1. Same labels as in Figure 3.

in very coarse case, they assume lower values in the range (10 < y + < 30) for PSS than LES. The wall-parallel RMS are globally underestimated.

Figure 5 shows the non-dimensional mean eddy-viscosity for LES and variance tensor components for PSS. The SGS eddy-viscosity and the wall-normal component a yy have similar profiles, and they are discussed together below.

They display common features for all the meshes used: in the viscous sublayer (y + < 5), they are practically null; in the buffer layer (5 < y + < 30), they rapidly increase and reach a peak in the range y + ∈ [10, 15], after which a smooth decay starts. In the log-law region (30 < y + < 150), the profiles for fine and coarse meshes decrease moderately and they eventually reduce to low values at the channel centre; the profile for very coarse mesh reports a more regular slope. For the coarse meshes, ν sgs shows slightly higher values than a yy near the channel center. This is possibly caused by the particular numerical implementation of the LES Samgorinsky model. However, this does not affect velocity statistics. Their values are moderate for the fine and the coarse grids (maximum 50% of the molecular viscosity), whereas they are of the same order of molecular viscosity for very coarse mesh. Hence, the SUS/SGS model plays a crucial role in this last case. The coordinates of the peaks correspond to the points where the damping starts. This is set to a fixed value y + B = 12.7 for the LU near-wall model, while it is variable for the van Driest model. With respect to the former, the latter is activated slightly closer to the wall in the fine mesh case, about at the same point in the coarse case, and slightly further from the wall in the very coarse case. In this last case, PSS provides a higher level of eddy-viscosity which reflects on the mean velocity and the streamwise RMS profiles (see Figures 3 and4), which are closer to the reference data in the buffer region. These results validate the pertinence of the LU wall-law model. to the fact that the wall-parallel component of variance tensor are not damped, but contribute to the energy dissipation term in equation ( 12). Dissipation is higher when the mesh degrades. The production terms are similar for PSS and LES: in the fine mesh case, they peak at y + 15 and y + 13 (respectively), in the coarse case they both peak at y + 19, while in the very coarse one at y + 35. It is worth to note that the PSS for the very coarse case yields a lower production close to the wall (5 < y + < 20), probably as an effect of the lower streamwise RMS (see also discussion Figure 5). The loss of energy due to SUS is only present in the pseudo-stochastic model; it assumes non-negligible negative values close to the wall (10 < y + ), and increases in magnitude as the mesh become coarser. It contributes to the total TKE dissipation. The turbulent compression/expansion term due to SUS is practically zero and does not contribute to the TKE budget.

Turbulent advection and compressibility

The additional terms that characterise the pseudo-stochastic model are here tc along wall distance. Average in time and wallparallel directions. Simulations with the three meshes described in Table 1.

< Φ + T C > y + < u + T A > y Ò Ó Ö× Ú ÖÝ Ó Ö×
strongly connected with the variance tensor behaviour in Figure 5. Figure 8 shows the Φ + tc isosurfaces of negative (blue) and positive (orange) near the bottom wall, at the last time configuration. They are organised in spots elongated in the streamwise direction, confined in the buffer and viscous layer. In accordance with the mean profile, positive spots are closer to the wall, while the negative ones are immediately above. The shape and the position of these structures suggest a possible correlation with the streaks turbulent structures (e.g. see [START_REF] Chernyshenko | The mechanism of streak formation in near-wall turbulence[END_REF]); however, an additional study is required to better investigate such a correlation.

Conclusion

The Finally, the pseudo-stochastic model is a promising alternative approach for turbulent modelling, that generalises the classical models and describes a richer physics. Mathematical and numerical investigations demonstrate the potential of the model.

  the profile is often extended till the half of the channel. Let us stress that these profiles are rigorously derived from the LU models. See[START_REF] Pinier | A model under location uncertainty for the mean velocity in wall bounded flows[END_REF] for an extensive validation on the pipe flow, turbulent boundary layer, and channel flows.

  ) is validated on turbulent plain channel flow at Re τ = 395. Subsequently, the pseudo-stochastic model is studied in detail on channel flow at Re τ = 590. Several simulations are performed changing the computational grid resolution, and the results of pseudo-stochastic simulation (PSS) are compared with a LES and the direct numerical simulation (DNS) of[START_REF] Moser | Direct numerical simulation of turbulent channel flow up to re τ =590[END_REF]. 4.1. Methodology and implementation Simulations are performed taking advantage of the open-source software OpenFOAM v6. This is a C++ library for computational fluid dynamics and uses the finite volume method. The LESs are carried out using the solver pisoFoam included in the standard software distribution. The implementation details can be found in the OpenFOAM documentation and in Jasak et al. (1999). The filtered classical Navier-Stokes equations are closed by the Smagorinsky model (24), with c s = 0.65. The van Driest function (25) for near-wall damping is used unless otherwise specified. The optimal parameters are set as n = m = 1, A + = 0.26, C δ = 0.158, which lead to a formulations similar to the original one by van Driest (1956). The PSSs are carried out using the home-made solver pseudoStochasticPisoFoam, developed by the authors at the Fluminance research group at INRIA Rennes (France). The pseudo-stochastic equations (4) are solved employing the Pressure-Implicit with Splitting of Operators (PISO) algorithm proposed by Issa et al.(1986) and[START_REF] Oliveira | An improved piso algorithm for the computation of bouyancy driven flows[END_REF]. The variance tensor is expressed by the isotropic constant model (13), corrected by the near-wall damping function (17) unless otherwise specified. The model constant is set to be c m = 2c 2 s in analogy with the Smagorinsky model. The damping parameter is set to be y + B = 2/κ = 12.7 after a theoretical estimation, confirmed by several test simulations. In order to regularise the damped profile of a yy , a smoothing filter is applied to the variance tensor.

  The fine meshes are such that the first cell is within y + = 1 and with 9 cells in y + ≤ 11, and the cell width in the wall-parallel plane are sufficient to ensure an accurate resolution of the boundary layer. The coarse and very coarse meshes still have a good vertical resolution but the streamwise and spanwise discretisation is reduced.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 the constant Smagorinsky SGS model, that has reached the statistical steady state.

Figure 1 :

 1 Figure 1: Non-dimensional mean velocity profiles along wall-normal direction for turbulent channel at Reτ = 395. Solid black, DNS by Moser et al. (1999). Top profiles: dash violet, analytical profile (17) derived from LU by Pinier et al. (2019). Bottom profiles: red symbols, PSS with near-wall model; red lines, LES with van Driest damping; blue symbols, PSS without near-wall model; blue lines, LES without van Driest damping.

Figure 2 :

 2 Figure2: Non-dimensional mean eddy-viscosity parameters along wall-normal direction for turbulent channel at Reτ = 395. Dash blue, LES with van Driest damping; solid red, PSS wall-parallel components of a ij ; circle red, PSS wall-normal component of a ij ; solid black, well-resolved LES in[START_REF] Armenio | A lagrangian mixed subgrid-scale model in generalized coordinates[END_REF] with spectral code.

  motion and minimising the model error along a fluid particle trajectory. The near-wall model has a crucial role in the correct damping of the eddy-viscosity close to the wall, both for LES and PSS. The LU near-wall model appears to accurately reproduce the slope in the region 5 < y + < y + B , while the van Driest model exhibits a larger deflection.

Figure 5 :

 5 Figure5: Non-dimensional mean eddy-viscosity parameters along wall-normal direction of turbulent channel at Reτ = 590. Simulations with the three meshes described in Table1. Solid line, sub-grid scale eddy viscosity from LES with van Driest damping; symbols, wallnormal variance tensor component from PSS with near-wall damping; dash lines, wall-parallel variance tensor components from PSS with near-wall damping.

Figure 6

 6 Figure 6 reports selected terms of the resolved TKE budget (12), averaged in time and wall-parallel planes, for the three meshes used. The time variation of TKE is made non-dimensional by the molecular viscosity. The equation for LES is obtained from (12) setting a ij = 0, except in the dissipation term where a ij = ν sgs δ ij in order to account for the dissipative effect of the sub-grid model. The dissipation profiles of the PSS and LES are identical except in the region y + < 20 close to the wall, where the PSSs have lower values. This is mainly due

Figure 6 :

 6 Figure6: Non-dimensional Turbulent Kinetic Energy (TKE) budget (12) along the wallnormal direction. Simulations with the three meshes described in Table1: from top to bottom: fine, coarse and very coarse mesh. Red lines with symbols, PSS. Blue lines without symbols, LES.

Figure 7 :

 7 Figure 7: Non-dimensional turbulent advection (wall-normal component) u +ta and nondimensional turbulent compressibility Φ + tc along wall distance. Average in time and wallparallel directions. Simulations with the three meshes described in Table1.

Figure 7

 7 Figure 7 displays the wall-normal component of non-dimensional turbulent advection u + ta,y = u ta,y /u τ and the non-dimensional turbulent compressibility Φ + tc = Φ tc ν/u 2 τ along wall distance. The other components of u ta are almost zero; thus, they are not reported. Globally, the magnitude of both quantities increases when the discretisation points decrease, since a larger part of the flow turbulence has to be modelled. In all the cases, wall-normal turbulent advection peaks at y + = 10 and is almost zero in the viscous sublayer and log-law region.The magnitude is quite small compared with the mean streamwise velocity: in the very coarse case, the peak of the vertical turbulent advection is 1.4% of the mean streamwise velocity at the same point. However, it generates a nonnegligible vertical velocity that drives the flow towards the wall. This qualifies u ta as a turbophoresis velocity, that advects the flow from a region of high to low turbulence level (quantified by the RMS velocity intensity). Turbulent compressibility presents a maximum at the end of the viscous sublayer y + = 5, and a minimum at y + B . It assumes moderate values. When positive (negative), it can be associated with a sort of fluid expansion (contraction) of the fluid due to turbulence.

Figure 8 :

 8 Figure 8: Positive and negative isosurfaces of Φ + tc near the bottom wall.Orange: isosurfaces Φ + tc = 0.5 max(Φ + tc ) = 0.013. Blue: isosurfaces Φ + tc = 0.5 min(Φ + tc ) = 0.0095.

  pseudo-stochastic model introduced byMémin (2014) is investigated both mathematically and numerically. Such a model is shown to be a generalization of the classical eddy-viscosity turbulent models, where the variance a ij plays the role of an eddy-viscosity tensor. Turbulence effects are not limited to energy dissipation, but induce additional phenomena as turbulent advection and compressibility that are not usually considered. Moreover, it does not rely on the restrictive physical assumptions related to Boussinesq's hypotheses. The turbulent kinetic energy budget is derived and presented, along with a nearwall model for a ij that is inferred from the analysis of boundary layer by[START_REF] Pinier | A model under location uncertainty for the mean velocity in wall bounded flows[END_REF]. A simple isotropic constant model for a ij is adopted for numerical simulation of turbulent channel flows, which are directly compared with an equivalent large-eddy simulation with the Smagorinsky model. In both cases, a near-wall damping function is used to correct the turbulent model in the proximity of the solid boundaries: the latter uses the classical van Driest function, the former employs the LU near-wall model here developed.First, the LU near-wall model is successfully validated in the channel flow at Re τ = 395. The eddy-viscosity tensor is correctly damped and exhibits a better agreement than the van Driest function with a reference solution obtained by a highly accurate simulations model. It is worth noticing that the LU model acts only on the wall-normal direction and depends on one single parameter (theoretically estimated); in contrast to the classical model that damps eddyviscosity for all the velocity components and requires to set several parameters.Hence, the former appears to impose a minimal correction and to have reduced empirical content. Second, the channel flow at Re τ = 590 is simulated using fine, coarse and very coarse meshes. Overall, the pseudo-stochastic simulations with the LU near-wall model are as accurate as the classical techniques. They show slightly better results when the computational grid is very coarse, since the van Driest model tends to excessively damp the eddy-viscosity. The PSS model is more effective in dissipating turbulent kinetic energy near the wall. A weak turbulent advection velocity is detectable between the viscous and buffer layer; such a velocity slightly advects the flow near the wall, form regions at high to low turbulent level (with respect to velocity RMS intensity). Hence, it is qualified as a turbophoresis phenomenon. In the same region, turbulent compressibility displays moderate positive and negative values; the visualisation of instantaneous isosurfaces suggests a possible link with the streaks turbulent structures.

Table 1 :

 1 Computational grid settings for numerical simulations of turbulent channel flow. The y + wall is the coordinate of the first point near the wall.

	Reτ Mesh	grid points	y + wall ÷ ∆y + max	∆x + ∆z +	λ
	395	fine	50 × 80 × 80	0.71 ÷ 25	50	23	5.00
		fine	96 × 96 × 96	0.71 ÷ 36	40	20	5.25
	590	coarse	64 × 64 × 64	1.14 ÷ 48	58	29	5.20
		very coarse 32 × 64 × 32	1.14 ÷ 48	116	58	5.20

Table 1 .

 1 Solid line, sub-grid scale eddy viscosity from LES with van Driest damping; symbols, wallnormal variance tensor component from PSS with near-wall damping; dash lines, wall-parallel variance tensor components from PSS with near-wall damping.

Table 1

 1

	0.02
	0
	0.02 -0.12 -0.1 -0.08 -0.06 -0.04 -0.02
	0.015
	0.01
	0.005
	0
	-0.005 -0.01 -0.015 0.1	1	10	100

: from top to bottom: fine, coarse and very coarse mesh. Red lines with symbols, PSS. Blue lines without symbols, LES.


Aymen Sendi
email: sendi@laas.fr

Grérory Besnard
email: besnard@laas.fr

Philippe Menini
email: menini@laas.fr

Chaabane Talhi
email: talhi@laas.fr

Frédéric Blanc
email: blanc@laas.fr

Bernard Franc
email: franc@laas.fr

Myrtil L Kahn
email: myrtil.kahn@lcc-toulouse.fr

Katia Fajerwerg
email: katia.fajerwerg@lcc-toulouse.fr

Pierre Fau
email: pierre.fau@lcc-toulouse.fr

Myrtil Kahan

Detection and Discrimination Methods of Sub-ppm Nitrogen Dioxide (NO2) with a Copper Oxide Sensor Operated with a Pulsed Temperature Modulation

Keywords: NO2, CuO nanoparticules, temperature modulation, gas sensor, selectivity

come   L'archive ouverte pluridisciplinaire

Introduction

Humans spend more than 90% of their time in closed environments, even though this indoor environment offers a wide variety of pollutants [START_REF] Lévesque | Indoor air quality[END_REF] [START_REF] Namiesnik | Indoor air quality (IAQ)". Pollutants, their sources, and concentration levels[END_REF]. Indoor air pollution is a real health threat, so measuring indoor air quality is important for protecting the health from chemical and gaseous contaminants. Nitrogen dioxide (NO2) is a dangerous pulmonary irritant [START_REF] Lévesque | Indoor air quality[END_REF]. NO2 is generated by multiple sources of combustion in indoor air, such as smoking and heaters, but it also comes from outside air (industrial sources, road traffic) [START_REF] Cadiergues | Indoor air quality[END_REF]. NO2 may have adverse effects of shortness of breath, asthma attacks and bronchial obstructions [START_REF] Koistinen | The INDEX project: Executive summary of a European Union project on indoor air pollutants[END_REF].

It is also classified as toxic by the "International Agency for Research on Cancer (IARC)" [START_REF] Loomis | The International Agency for Research on Cancer (IARC). Evaluation of the carcinogenicity of outdoor air pollution: focus on China[END_REF], hence the necessity for sensor development for accurate NO2 detection is an acute need. Among sensors techniques, the metal oxide gas (MOX) sensors are promising candidates because of their high performance in terms of sensitivity on top of their low production cost. The copper oxide (CuO) material is highly studied because of its high sensitivity and its ability to detect oxidant gaseous compounds, but also for other indoor air pollutants, such as acetaldehyde (C2H4O), formaldehyde (CH2O), NO2, CO, etc. However, CuO (as well as other metal oxides) suffers from a major disadvantage which is the lack of selectivity with respect to targeted gas.

In this study, our main objective is to develop an innovative and simple pulsed-temperature operating mode associated with an efficient data processing technique, which enables good selectivity toward NO2 in gas mixtures (combination of two gases). This technique is based on few parameters extracted from the dynamic response of sensor versus temperature changes in a gaseous environment. These parameters are: the normalized sensing resistance, the initial slope, the intermediate slope and the final slope of the response of NO2 against different reference gases, such as C2H4O, CH2O and moist air.

These four parameters allow us to follow the synthetic evolution of the adsorption phenomenon of gaseous molecules in the surface of the metal oxide during the gas injection. This synthetic evolution is different for each gas which allows us to define detection meanings for each gas in the range of the defined concentration. This technique greatly facilitates the discrimination of NO2 toward other referent gas and improves the selectivity of CuO under different gas contaminant and condition in indoor air.

The selectivity of NO2 was examined in relation to air moisture with 30% humidity, C2H4O at 500-ppb, CH2O at a concentration of 500-ppb and the binary mixture of these gases with 300-ppb of each.

In Section II of the paper, we describe the materials and methods used in our work. The Section III presents our results and the discussion. We conclude this work in Section IV.

Materials and Methods

Description of the Metal Oxide: the Copper Oxide

CuO nanoparticles have been obtained by the reaction of ambient air during 24 hours on copper acetamidinate precursor (Cuamd) in the presence of 5 molar equivalents of octylamine (OA). [START_REF] Jońca | Organometallic synthesis of CuO nanoparticles: Application in low temperature CO detection[END_REF] Excess of OA in the medium is removed by several precipitation stages by centrifugation and washing with THF. In Figure 1 we show a SEM micrographic of a CuO. The sensitive layer made of CuO nanoparticles is deposited by inkjet [START_REF] Morozava | Thin TiO2 films prepared inkjet printing of the reverse micelles sol-gel composition[END_REF] printing on a silicon microhotplate [START_REF] Dufour | Technological improvements of a metal oxide gas multi-sensor based on a micro-hotplate structure and inkjet deposition for an automotive air quality sensor application[END_REF] presented with more detail in the next section. We chose this technique because it is well adapted for various material integration on multisensor platform (our final goal) without additional photolithographies Moreover, it has been less studied toward the other deposit technique [START_REF] Kukkola | Novel printed nanostructured gas sensors[END_REF]. The ink is prepared with 5% CuO weight, which was dispersed in ethylene glycol by an ultrasonic bath for about one hour. The dispersions obtained were allowed to settle for 24 hours. The final ink was collected and then used for printing using Altadrop equipment, where the numbers of the deposited drops of ink were controlled [START_REF] Dufour | Technological improvements of a metal oxide gas multi-sensor based on a micro-hotplate structure and inkjet deposition for an automotive air quality sensor application[END_REF]. The principle of operation of this device and the control of the number of drop is described in Figure 2. This technique is simple and allowed us to obtain reproducible thicknesses from one to a few micrometers depending on the number of deposited drops. In addition, this technique permits to have a precisely localized deposit without need of additional complex photolithographic steps. The CuO layer is finally deposited on the microhotplate and annealed in ambient air from room temperature to 500°C (rate 1°C/min) followed by a plateau at 500°C for 1 hour before cooling to room temperature (1°C/min). This initial temperature treatment is necessary because CuO requires operating temperatures between 100°C<T<500°C. The thermal pretreatment is necessary to generate ionized oxygen species in atomic or molecular form at the oxide surface and therefore to improve the reactivity between the reacting gas and the sensor surface [START_REF] Jońca | Organometallic synthesis of CuO nanoparticles: Application in low temperature CO detection[END_REF].

Microhotplate

The microhotplate has been designed in the technology center of LAAS-CNRS in Toulouse in France, and optimized to increase the temperature up to 600°C with low energy consumption (only 60 mW), a good reliability and high robustness [START_REF] Yoboue | An Industrializable silicon based microhotplate for metal oxide gas sensors[END_REF].

This microhotplate is designed on a silicon mass present a support for a bilayer membrane composed by SiO2/SiNx on which the heater and the sensitive electrode are placed. The metallization of these electrodes is composed by a bilayer Ti/Pt then insulated by a SiO2 passivation layer. In Figure 3 we see a top view and a cross sectional view with different physical layer.

Operating Mode

In this study, we have used a pulsed temperature profile, which showed that optimized sensitivity can be achieved with the use of two different temperature stages at 100°C and 500°C respectively. This dual temperature protocol also reduces the total power consumption of the device. (See Figure 4) In all these tests we applied these two forms of modulation in the same profile in order to see the effect of the transition from a high temperature to a low temperature on the CuO sensitivity, passing through intermediary temperature without return to the step of the high temperature (in the case of pulsed down temperature modulation) and with a return towards this level in order to see the effect of heating the CuO follows a low temperature step on the gas detection. In Figure 6 we show the results of these four tests.

Each bar shows the results of the normalized resistance of a gas injection for each time step by applying on the heater the temperature profile presented above. By comparing the sensitivity bars and the transient responses during each temperature step, we concluded that the best performances under NO2 was obtained with the pulsed down temperature modulation form and with a 30 seconds duration step. The transition that carry out the best sensitivity is from the highest temperature (500°C) to the lower temperature (100°C) (grey color bar in fig6), that explains our choice for the final temperature profile applied on our CuO sensor. CuO needs a high temperature at 500 ° C for cleaning the surface of the metal oxide and a low temperature at 100°C for measuring the adsorption under gas. In this way we obtained the optimal temperature profile for the CuO presented in figure 4.

This profile has then be applied on the heater of our CuO sensor for all further tests.

Material and Methodology

The CuO sensor was placed in a 250 ml test chamber and the test conditions were as follow: -A flow rate of 200 ml/min, controlled by a digital flowmeter.

-A relative humidity (RH) level of 30% is obtained by bubbling synthetic air flow controlled by a mass flow controller.

-The measuring chamber is at ambient temperature, controlled by a temperature sensor placed inside the vessel.

-A bias current is applied to the sensitive layer, controlled by a Source Measure Unit (SMU). We started a test with single gas injections after a phase of two-hours stabilization in humid air, then with injections of binary mixtures for 16-min. During 32-min, moist air is injected between two successive gas injections. This time is enough to clean the chamber and stabilize the sensor to its baseline.

The gas injections concentrations are summarized in Table 1.

A schematic representation of these injections is presented in Figure 7.

Table 1. Injected gases concentrations [START_REF] Sendi | Sub-ppm Nitrogen Dioxide (NO2) Sensor Based on Inkjet Printed CuO on Microhotplate with a Pulsed Temperature Modulation[END_REF]. During all the experience (stabilization phase, gases injections and stage between two successive gases) , the sensor is powered by a square signal voltage applied on the heater in order to obtain two temperature steps as shown in Figure 4. To ensure a constant overall flow, we adapted the gas injection sequences duration in correlation with the heating signal period (see Figure 8).

Gas

The resistance variation of CuO is measured under a fixed supply current of 100 nA in order to obtain voltage measurements in the range of few volts, far from compliance limit of the measuring device (20V). We also verified that this procedure (temperature cycling) doesn't affect the sensor reproducibility in terms of baseline or the sensor sensitivity.

Under such test conditions, we achieved a continuous 6.5 hours testing period without observing any drift on the raw sensor response. The sampling period is 500 ms, which gives us 60 points on a 30second response step, this acquisition rate being enough for accurate data processing [11]. Fig. 8. Diagram of a gas sequence [START_REF] Sendi | Sub-ppm Nitrogen Dioxide (NO2) Sensor Based on Inkjet Printed CuO on Microhotplate with a Pulsed Temperature Modulation[END_REF].

Finally, we analyzed the sensor responses on each steps, according to the different gas injections, using a simple method of data processing in order to have the better selectivity of NO2.

Results and Discussion

Method of Analysis

During each gas injection, 16 periods of temperature modulation are applied. After verifying the reproducibility of sensor responses along these cycles, we only present here the responses of the last cycle, which is stabilized and reproducible from one cycle to another. As mentioned previously, we used new simple data treatment methods to obtain the better selectivity toward NO2 with respect to several interferent gases.

Among the multiple possible criteria, we chose representative variables that take into account the dynamic sensor behavior during a change of gaseous conditions and during a pulsed temperature; these criteria are obtained from the sensor resistance slopes during the gas response on each 30-second-steps. The data acquisition relies on the decomposition of the response into three distinct domains in order to follow the synthetic evolution of CuO toward the different gas injections. (See Figure 9) These domains are: -Starting Slope: from the 1 st point to the 10 th point (in yellow).

-Intermediate slope: from the 10 th point to the 30 th point (in red).

-Final slope: from the 30 th point to the 60 th point (in black).

The normalized resistance is measured from the last point on each step: this is the absolute difference between the resistance of the sensor under a reference gas (like moist air) and the resistance of the sensor under targeted gas(es), at the final cycle of each injection: By treating these four parameters, a good selectivity of NO2 can be obtained with respect to moist air, C2H4O and CH2O.

𝑅𝑅𝑅𝑅 = � 𝑅𝑅𝑅𝑅𝑅𝑅𝑅𝑅-𝑅𝑅𝑅𝑅𝑅𝑅𝑅𝑅 𝑅𝑅𝑅𝑅𝑅𝑅𝑅𝑅 � * 100 (1)

Results

Before starting the data treatment from each slope, we have based on the raw data and the graph of the gas injection responses. Figure 10 shows the graph of the last cycle of a gas injection with the four temperature steps present the temperature profile applied. It is clear from this graph that the initial, the intermediate and the final slope at 500 ° C are flat and near to zero. (See Figure 10)

Fig.10. Example of a raw response gas injection

In the next we treat the three domains one by one

Initial Slope

The slope at the origin is calculated on the first 10 points of each temperature step of the 8 th cycle of each injection. Regarding the reference gas (moist air); we took the response of the 8th cycle of the last sequence under humid air before the injection of gases. The values of these slopes (in Ohms/ms) are shown in Figure 11; each bar represents the value of the slope at the origin of each gas at 500 and 100°C. Fig. 11. Representation of the slopes at the origin under different gases of the 8 th cycle in "Ohms/ms" according to the temperature of the sensors [START_REF] Sendi | Sub-ppm Nitrogen Dioxide (NO2) Sensor Based on Inkjet Printed CuO on Microhotplate with a Pulsed Temperature Modulation[END_REF].

The Figure 11, clearly shows that the calculations of this parameter enable us to differentiate NO2 from the other reference gases by measures on the step at 100°C. Regarding the other step at 500°C we note that the response is almost zero, because the transition from a cold to a hot state decreases the detection sensitivity and therefore reduces the sensors resistance variations under gas. We also note that with this criterion we evidence a significant difference between the value of the starting slope under NO2 compare to others under other gases and mixtures without NO2 [START_REF] Sendi | Sub-ppm Nitrogen Dioxide (NO2) Sensor Based on Inkjet Printed CuO on Microhotplate with a Pulsed Temperature Modulation[END_REF].

Intermediate Slope

The intermediate slope is calculated from the 10 th point to the 30 th point of each temperature step of the last cycle during each gas injection and compared with the value under air. The calculated values are shown in Figure 12. According to Figure 12, the value of the intermediate slope of individual NO2 injection or of the gas mixtures in which there is NO2, is not prominent compared with the other reference gases even on the plateau at 100°C. This parameter is less effective than the previous one to detect NO2 in gas mixtures [11].

Final Slope

The final slope is the slope of the second half of the gas response; it is calculated between the 30 th point and the 60 th point. The response to the different temperature steps of the last cycle during NO2 injection and the reference gases is presented in Figure 13 [11].

Fig.13. Representation of the final slopes under different

gases at the 8 th cycle in "Ohms/ms" according to the temperature of the sensors [START_REF] Sendi | Sub-ppm Nitrogen Dioxide (NO2) Sensor Based on Inkjet Printed CuO on Microhotplate with a Pulsed Temperature Modulation[END_REF].

It is worth noting that the value of the final slopes of NO2 or gaseous mixtures, which contain the NO2, are very different compared to the other reference gases on the two stages at 100°C only. This parameter allows us to select NO2 with respect to other interfering gases.

Normalized Resistance

As previously presented, the normalized resistance is calculated with respect to humid air from each last cycle level. The reference resistance used is the resistance of each stage of the last wet air sequence before the gas injection. The results obtained from these calculations are shown in Figure 14. The Figure 14 shows a slight variation in the values of the normalized resistance between two similar and successive temperatures for the same gases. This slight variation can be explained by data dispersion which is +/-2%, due to the fact that the normalized resistance is calculated from the raw values during the gas injection and the raw values of resistances during the injections of humid air that may be slightly different between two similar and successive temperature stages.

The CuO sensor response to sub-ppm NO2 levels when injected individually or in combination with another gas, is specific when measurements are taken on the low temperature plateau at 100°C [11].

Conclusions

The selectivity and sensitivity of our CuO sensor has been studied by different operating modes and simple methods of analysis. Specific temperature modulation was applied to the metal oxide with the use of temperature steps at 500 and 100°C. The response of CuO sensitive layer toward gases representatives of indoor air pollution (C2H4O, CH2O, NO2, humid air) has been studied. These responses were analyzed with several parameters, such as the study of the slope of resistance variation at the origin, the intermediate slope, the final slope and the normalized resistance measured at each temperature steps. By using these 4 parameters we can to define a new PCA, different from what exists in the literature of gas sensors (PCA with only the study of standardized resistance). The effect of adding new elements to the raw data processing of a response improves the selectivity of the detected gases.

The study of these different parameters shows that the CuO material is able to detect sub-ppm levels of NO2 with a good selectivity compared to different interfering gases. This study shows new parameters applicable to the PCA on any metal oxide in order to improve its selectivity.

To still improve the selectivity of gas sensor device to a larger variety of polluting gases, we plan to integrate these CuO sensors in a multichip system, which will allow us to use in parallel new metal oxide layers with specific temperature profiles and data analysis criteria.

Fig. 1 .

 1 Fig.1. SEM micrographic of a CuO.

Fig. 2 .

 2 Fig.2. Schematic diagram of the inkjet machine.

Fig. 3 .

 3 Fig.3. Micro-hotplate gas sensor: (a) Top view, (b) Cross sectional view.

Fig. 4 .Fig. 5 .

 45 Fig.4. CuO temperature profile.This optimal profile was obtained thanks to a specific study on the during of each step and a different temperature modulation forms in order to obtain the best sensitivity of the copper oxide under NO2.In this study we launched four tests by changing the duration of the step for each test and injected a 300ppm of NO2 to see where we can obtain the better sensitivity:1-1 th test : with time step equal 2 seconds 2-2 sd test : with time step equal 10 seconds 3-3 th test : with time step equal 30 seconds 4-4 th test : with time step equal 60 seconds For each test we applied a temperature profile with two modulation forms: -Pulsed down temperature modulation with a return to 500 ° C after each pulse by applying on the heater successively the following temperature: 500 -400 -500 -300 -500 -200 and 500 -100°C. -Descending step temperature modulation by applying on the heater successively the following temperature: 500 -400 -300 -200 and 100 ° C. (See Figure5)

Fig. 6 .

 6 Fig.6. Representation of the normalized resistance of CuO sensor under 300ppb of NO2 at different temperature stages.

Fig. 7 .

 7 Fig.7. Synoptic representative of a sequence of gas injections [11].

Fig. 9 .

 9 Fig.9. Representative diagram of a response of a gas sensor during an injection cycle and showing the 3 domain slopes [11].

Fig. 12 .

 12 Fig.12. Representation of the intermediate slopes under different gases at the 8 th cycle in "Ohms/ms" according to the temperature of the sensors [11].

Fig. 14 .

 14 Fig.14. The normalized resistance of the different gas injections at the 8 th cycle according to the temperature [11].

Sensors & Transducers, Vol. 0, Issue 0, Month 2016, pp.

Acknowledgments

The authors express their gratitude to neOCampus, the university project of the Paul Sabatier University of Toulouse for the financial support and the Chemical Coordination Laboratory of Toulouse for the preparation of the CuO nanoparticle powder. This work was also partly supported by the french RENATECH network.