Wassim Zeidan

Nicolas Mazellier
email: nicolas.mazellier@univ-orleans.fr

Emmanuel Guilmineau

Pierre-Yves Passaggia

Azeddine Kourta

P.-Y Passaggia

A near-wake survey of an Ahmed body comparing low-& high-fidelity numerical models with experiments

Keywords:

A full-scale square-back Ahmed body is analysed experimentally and numerically for a freestream flow velocity of 45m/s corresponding to a Reynolds number = 7.7×10 5 . Numerical simulations are performed for matching inflow conditions and different turbulence models. The main goal of this study is to examine the effectiveness of the Delayed Detached Eddy Simulations and the Reynolds Averaged Navier Stokes -SST turbulence model to capture the mean flow physics. Local and global aerodynamic parameters are compared with experimental results, including boundary-layer parameters evolving over the roof, the wake characteristics, the shear layer, and the mass entrained into the wake. A particularly good agreement is found between the experiment and the DDES. An in-depth analysis of the shear-layer growth, mass-entrainment and momentum budget analyses along the recirculation region interface, allow for identifying and quantifying the differences in predicting drag and the mechanisms at stake between the different approaches. The momentum budget, performed on the experimental data shows that pressure can be accurately obtained from experiments and closely follows that computed from high-fidelity numerical simulations. These results confirm the effectiveness of DDES in studying flows at high Reynolds numbers. The mass/momentum budget along the recirculation-region interface is also found as a possible mean to quantify road vehicle's aerodynamic performances, naturally extending to more complex geometries.

Introduction

The transportation industry, and in particular, road vehicles, remain one of the main contributors to global energy consumption and pollutant emissions. Car manufacturers must therefore comply with environmental regulations such as CO 2 -emissions standards which are becoming increasingly strict. Reducing the aerodynamic drag through flow control strategies stands as an appealing solution to comply with these policies but performing statistically converged numerical simulations at such high Reynolds numbers remains a challenge. [START_REF] Ahmed | Some salient features of the time-averaged ground vehicle wake[END_REF] introduced a simplified yet ubiquitous vehicle design, akin to a wide range of commercial road vehicles shapes. The objective of various experimental works including [START_REF] Spohn | Flow separations generated by a simplified geometry of an automotive vehicle[END_REF]; [START_REF] Beaudoin | Cavitation as a complementary tool for automotive aerodynamics[END_REF]; Conan, Anthoine and Philippe (2011); Grandemange, Gohlke and Cadot (2013c); [START_REF] Volpe | Experimental characterization of the unsteady natural wake of the full-scale square back ahmed body: flow bi-stability and spectral analysis[END_REF], numerical studies such as [START_REF] Krajnović | Influence of floor motions in wind tunnels on the aerodynamics of road vehicles[END_REF]; Lucas, Cadot, Herbert, Parpais and Delery (2017); [START_REF] Hesse | Simulation of wake bimodality behind squareback bluff-bodies using les[END_REF] regarding square-back Ahmed bodies, and Serre, Minguez, Pasquetti, [START_REF] Serre | On simulating the turbulent flow around the ahmed body: A french-german collaborative evaluation of les and des[END_REF]; [START_REF] Guilmineau | Assessment of turbulence models for flow simulation around the ahmed body[END_REF] for slanted Ahmed body is the study of the near-wake flow at high Reynolds number. Nevertheless, reattachment happens relatively fast since the frontal body shape is designed to avoid a large separation [START_REF] Ahmed | Some salient features of the time-averaged ground vehicle wake[END_REF]. The afterbody separation generates a very large recirculation bubble downstream of the body, contributing to the dominant part of the drag, knowing that 50% to 70% of the total drag is produced by the rear part of the vehicle. From this issue, we focus our analysis on the near-wake region behind an Ahmed body. The aim is to provide a generic approach to diagnose drag and address momentum budgets in the case of three-dimensional bluff bodies.

The near wake is characterised by a recirculation region which can be delimited from the outer flow using two different flow characteristics; ( ) the recirculating streamline closest to the outer flow, known as the Dividing Streamline Interface (DSI) [START_REF] Passaggia | Aerodynamic drag modification induced by free-stream turbulence effects on a simplified road vehicle[END_REF] shown in black in figure 1 and( ), the surfaces where the streamwise velocity is null [START_REF] Stella | Scaling of separated shear layers: an investigation of mass entrainment[END_REF], known as the Recirculating-Region Interface (RRI) shown in red-blue in figure 1. Note that in the case of a wake, both the DSI and the RRI originate and terminate at the separation and the saddle point closing the recirculation bubble respectively. These generic interfaces can therefore be used to perform momentum budgets and link the topology of the near wake and the physical quantities of the flow to the forces acting upon the vehicle. Characterising drag can be therefore linked simultaneously to a geometric problem (i.e. the topology of the near wake) and a physical problem (i.e. estimating mechanical loads). The geometric problem aims at understanding the topological features of the wake. A particularly important quantity for a momentum budget analysis is the local slope of the interface which, in two dimensions, links the spatial derivatives and can be readily obtained from differential calculus as d = d ∕d

| | | d = d
, where is the local interface's slope and is the surface of the RRI or the DSI. The slope allows for linking physical quantities projected along the vertical direction with the horizontal contributions in the budget and reduces the momentum budget to a single equation linking the topology to the physical quantities. An interesting feature of both the RRI and the DSI is that the mean flow contributions naturally cancel when considering the budget for the drag force acting at the base of the vehicle. Whence, the budget along either the DSI and the RRI in the case of a 2D flow reduces to

, = ∫ + ⟨ ′2 ⟩ d + ∫ ⟨ ′ ′ ⟩ d = ∫ + ⟨ ′2 ⟩ + ⟨ ′ ′ ⟩ d ,
where , is the pressure contribution to the drag force at the base of the vehicle per spanwise unit length, is the mean pressure, while primed quantities denote turbulent quantities. Note that the RRI allows for computing the mass entrained across the recirculation region while this quantity naturally cancels when considering the DSI. In addition, the RRI is well defined for three-dimensional flows, contrary to the DSI. In this paper, we therefore use the RRI in order to analyse the flow in the near wake of an Ahmed body at a high Reynolds number for both numerical simulations and laboratory experiments for a quantitative comparison.

A remarkable body of work regarding high-fidelity simulations of the flow around squared-backed Ahmed body at high Reynolds number has been reported by [START_REF] Lucas | A numerical investigation of the asymmetric wake mode of a squareback ahmed body effect of a base cavity[END_REF]; Dalla Longa, [START_REF] Dalla Longa | Simulations of the bi-modal wake past three-dimensional blunt bluff bodies[END_REF]; [START_REF] Evstafyeva | Simulation and feedback control of the ahmed body flow exhibiting symmetry breaking behaviour[END_REF]; [START_REF] Hesse | Simulation of wake bimodality behind squareback bluff-bodies using les[END_REF], regarding large-eddy simulations and wall-modelled large-eddy simulations. In these studies, the authors analysed the role of the modelling strategy on the the dynamics. Nevertheless a thorough analysis comparing experiments and simulations of the same body is still lacking. In particular, systematic approaches including global budgets to analyse the near-wake features have yet to make consensus. Thus, the importance of the simulation presented in this article is the high Reynolds number, close to one million, which replicates a wind-tunnel experiment characteristics from the incoming boundary layer to the geometry of the near wake. At such high Reynolds numbers, well-resolved large-eddy simulations [START_REF] Dalla Longa | Simulations of the bi-modal wake past three-dimensional blunt bluff bodies[END_REF] are out of scope since the necessary simulation time to captures several wake reversals is too computationally intensive. Similar to [START_REF] Hesse | Simulation of wake bimodality behind squareback bluff-bodies using les[END_REF], we use an intermediate approach such as delayed-detached-eddy simulations [START_REF] Guilmineau | Assessment of turbulence models for flow simulation around the ahmed body[END_REF][START_REF] Guilmineau | Assessment of hybrid rans-les formulations for flow simulation around the ahmed body[END_REF] over extensive simulation times to make sure that the temporal statistics of the mean flow are well converged. Besides, coupling between experimental and simulation analysis should allow for a more precise understanding of the flow physics. Thus, we report a complete analysis of the boundary layer developing over the roof, together with a mass and a momentum budget in the near wake in order to evaluate the quality of the mean field between experiments, high-fidelity numerical simulations and a standard Reynolds-averaged Navier-Stokes approach, using the RRI for the analysis.

The paper is organised as follows, §2 provides details about the experimental setup, §3 details the numerical strategy, and §4 provides a summary of the mathematical strategy to analyse the flow. Results are summarised and compared in §5 while conclusions are drawn in §6.

Experimental Setup

The experiment was performed at PRISME laboratory in Orleans University (France) in the wind tunnel of test section 2×2 m 2 wide, and 5 m long. The model is a full-scale Ahmed body [START_REF] Ahmed | Some salient features of the time-averaged ground vehicle wake[END_REF], with length = 1044 mm, width = 389 mm, height = 288 mm, ground clearance ℎ = 45 mm (15% of the model height), within a blockage of 2.8%. The model was mounted via a plate located 0.47 m from the ground. The -axis is aligned with the streamwise direction, the -axis is vertical and positive upwards (as shown in figure 2). The origin is at the bottom of the base of the model at the mid wheelbase. An external strain gauge balance was used to measure the aerodynamic forces, located under the working section of the tunnel. The chosen sampling frequency is 1kHz but the dynamic response of the balance is 1Hz. The Ahmed body is linked to a six-components strain gauge balance through its feet, avoiding any contact with the ground. We focus here on the drag force . The drag coefficient is measured with a 0.5% precision and is well within the range reported in the literature for a similar geometry and flow condition [START_REF] Volpe | Experimental characterization of the unsteady natural wake of the full-scale square back ahmed body: flow bi-stability and spectral analysis[END_REF][START_REF] Lucas | A numerical investigation of the asymmetric wake mode of a squareback ahmed body effect of a base cavity[END_REF]. Note that the same setup was studied in [START_REF] Passaggia | Aerodynamic drag modification induced by free-stream turbulence effects on a simplified road vehicle[END_REF] to analyse the effect of free-stream turbulence on the wake but for a slightly lower free-stream velocity and details about the experimental setup are shown in figure 1.

The pressure at the base was measured by covering the back-square face with a grid of 81 pressure taps connected through 600 mm long vinyl tubes, with an internal diameter = 1.56 mm, and with a manufacturer quoted accuracy of ±6 . Furthermore, pressure taps were located with a fine distribution close to the surface edges. In addition, two MicroDAQ pressure scanners were used with a sampling frequency of 50Hz. Note that the temperature in the test section was 38 C and that the kinematic viscosity of air was computed accordingly. In addition, two dimensional PIV was performed using a LaVision system to analyse the topology of the near wake in themid plane, at =0, the image size is 392 × 603 mm 2 with a spatial resolution of 2 mm for each vector. Measurements were acquired at a frequency of 2.49 Hz. The data were processed using the DPIVSoft toolbox (see [START_REF] Passaggia | Estimating pressure and internal-wave flux from laboratory experiments in focusing internal waves[END_REF] and [START_REF] Meunier | Analysis and treatment of errors due to high velocity gradients in particle image velocimetry[END_REF]). Further details regarding PIV can be found in [START_REF] Passaggia | Aerodynamic drag modification induced by free-stream turbulence effects on a simplified road vehicle[END_REF].

Numerical simulations

Flow solver

The ISIS-CFD flow solver, uses an incompressible unsteady Reynolds-averaged Navier-Stokes (URANS) method. The solver is based on a finite-volume method to build the spatial discretization. The unstructured discretization is face based with adaptive mesh refinement strategies. A second-order backward difference scheme is used for the temporal discretisation. A centred scheme is used for the diffusion terms whereas the convective fluxes use a blended scheme (i.e. 95% central differences and 5% upwind differences for the hybrid RANS-LES model, or 80% central and 20% upwind for the RANS model [START_REF] Krajnović | LES of the flow around a generic wheel in a wheelhouse[END_REF]). Two fully unstructured hexahedral meshes were constructed using HEXPRESS, the first with 32 million cells and the second with 56 million cells. The convergence was chosen referring to physical parameters such as the drag coefficient, the recirculating region length alongside the experiment.

To compromise between the computational costs and accuracy, a mesh study was conducted and the first mesh was considered. Snapshots are shown in figure 3b for illustration. In this paper, two turbulence models are used: the RANS k-SST [START_REF] Menter | Two-equation eddy-viscosity turbulence models for engineering applications[END_REF] and the Delayed Detached Eddy Simulation (DDES). The DDES model is a SST-based DDES proposed by [START_REF] Gritskevich | Development of DDES and IDDES Formulations for the k-Shear Stress Transport Model[END_REF]. The benefits of DDES is the ability to split the computational domain into two zones; Near-wall regions are covered by the first zone where the standard RANS equations are resolved while a LES approach is used in the second zone. The velocity field is obtained from the momentum conservation equations and the pressure field is extracted from the mass conservation constraint or continuity equation, transformed into a pressure-equation. The pressure equation is obtained by the [START_REF] Rhie | A Numerical Study of the Turbulent Flow Past an Airfoil with Trailing Edge Separation[END_REF] interpolation. The momentum and pressure equations are solved in a segregated manner as in the simple coupling procedure [START_REF] Issa | Solution of the implicitly discretised fluid flow equations by operator-splitting[END_REF].

The same simplified vehicle model (same as used in the experiments) is discretised on a computational domain with dimensions 13.2 × 1.7 × 1.9 (see figure 3a). The upstream distance between the inlet and the body were chosen, and at a position 1.5 upstream the front part of the body. A boundary layer measurements is taken and compared with experiment. Figure 4 shows the boundary-layer comparison between experiments and simulation at a distance 1.5 in front of the Ahmed body. The figure shows a comparable result between the simulation and experiment at velocity equal to 45 m/s, which validates the size of the computational domain.

Boundary conditions were chosen as follows, a 45 m/s inlet was imposed, a no-slip condition was considered for the body and the road to enable boundary-layer growth. In order to capture near-wall features, a refinement was performed defining three boxes covering the whole body and zooming on the rear part in which the cell length was 1.8mm. The coordinate axis system was defined similar to the experiment to facilitate the comparison (see figure 3a), with an origin at the mid-base bottom part of the body. The non-dimensional time step Δ ⋆ = 0.1 captures the boundary layer correctly and is higher than the time step calculated by Δ ⋆ = Δ ∞ ∕ , with Δ = Δ ∕ , where is the average velocity in the wake region. Note that 480 processors were used, with a computational time ⋆ = 1250, where ⋆ = ∞ ∕ , is the physical time in seconds, ∞ is the free-stream velocity and is the height of the body. The resulting mean flow is shown in figure 3c while a snapshot of the turbulent flow is depicted in figure 3d by means of the -criterion coloured with the instantaneous horizontal velocity.

Mathematical tools

In what follows, the flow quantities are non-dimensionalized using the height of the body and the free-stream velocity ∞ such that

⋆ = ∞ ∕ , ⋆ = ∕ , ⋆ = ∕ ∞ , and ⋆ = ∕ 2 ∞ , (1)
unless specified otherwise. Note that is the density of the fluid. Further, the analysis follows a Reynolds-Averaged Navier-Stokes (RANS) decomposition for the analysis between a mean and a fluctuation such that

( , ) = ( ) + ′ ( , ), ( ) + ′ ( , ) .
(2)

The time-averaged Navier-Stokes that is used throughout this study writes in non-dimensional form where = ∞ ∕ is the Reynolds number kept nearly constant at ≈ 7.7 × 10 5 ± 1% throughout this study. Note that the viscous term on the right-hand side of equation ( 3) is negligible and we do not consider its contribution hereinafter.

Mean pressure reconstruction of ( , ) is performed in the experiment in order to quantify the role of the mean pressure and Reynolds stresses in a momentum budget. The procedure first solves (3), which is akin to solving the Poisson problem

∇ 2 ⋆ = -⋆ ⋆ ⋆ ⋆ + ⟨ ′ ′ ⟩ ⋆ , (4)
but in a least-square sense [START_REF] Harker | Least squares surface reconstruction from measured gradient fields[END_REF][START_REF] Harker | Direct regularized surface reconstruction from gradients for industrial photometric stereo[END_REF][START_REF] Passaggia | Estimating pressure and internal-wave flux from laboratory experiments in focusing internal waves[END_REF], to obtain a first guess of the static pressure 0 and then uses the pressure taps at the base of the body (at = 0) to compute the static pressure. The true pressure ⋆ is recovered using the adjoint method derived in Shanmughan, Passaggia, Mazellier and Kourta (2020), and adapted from [START_REF] Passaggia | Aerodynamic drag modification induced by free-stream turbulence effects on a simplified road vehicle[END_REF] for the boundary conditions. Note that in the case of a symmetrised flow with respect to = 0, the missing out-of-plane components of the velocity for the PIV are not problematic since the mean transverse ′ -component of the flow is null by construction in the ( -) plane for the symmetrised flow while the transverse Reynolds-stress tensor gradient (i.e. ⟨ ′ ′ ⟩ ⋆ ) has a negligible contribution in the recirculation region as will be shown later. The pressure reconstruction therefore makes little assumptions and is essentially dependent on the quality of the measurements.

Results and comparisons

Boundary-layer analysis

The boundary layer thickness, can be defined as the distance from the wall to the point where the mean streamwise velocity approaches 99% of the local free-stream velocity ∞ . More interesting for this comparison is the displacement and the momentum thicknesses defined as

⋆ = ∫ ∞ 0 1 - ∞ and = ∫ ∞ 0 ∞ 1 - ∞ , ( 5
)
where ⋆ is linked to the flow-rate deficit inside the boundary layer while is directly related to the skin-friction coefficient = 2 ∕ 2 ∞ . The shape factor is given by = ⋆ ∕ and allows for a first assessment of the state of the boundary layer. The Reynolds number based on the momentum thickness = ∞ ∕ while the non-dimensional parameter for the pressure gradient writes = ( ⋆ ∕ ) ⋆ ∕ ⋆ . The displacement and momentum thickness were evaluated using a trapezoidal rule. The boundary layer measurements were specified 12 cm (i.e. ⋆ = -0.7) upstream the top trailing edge of the body, over the roof. The origin is set at the mid-bottom of the body base, shown in figure 1. Hence, several points were acquired to measure the boundary-layer properties using a single-component hot wire with a sampling frequency of 60 kHz. However, due to the limitation in the experimental setup, a minimum distance of 0.5 mm from the wall remained insufficient to enter the viscous sub-layer. Experimental and simulation results are reported in table 1. The boundary layer can be considered as fully turbulent, since both conditions > 3000 and > 1.2 are met but it should be noted that the boundary layer may have not reached equilibrium since ≈ 1.3. Figure 5a shows the mean streamwise velocity as a function of the wall-normal distance from the roof of the body normalized with the free-stream velocity and the boundary layer thickness respectively. A good match between experiments and the simulations (for both turbulence models) is observed. A slight difference is noticed along the evolution of the streamwise velocity profile between the experiment and simulations (see figure 5a). In a similar study, [START_REF] Trip | Bluff body boundary-layer modification and its effect on the near-wake topology[END_REF] report velocity profiles distinctions because of the state of the turbulent boundary layer which is expected since ≈ 1.3. As a result, both cases (experiments and simulation) demonstrate a turbulent boundary layer.

The inner region of the boundary layer is defined from = 0 up to = 0.1 [START_REF] Pope | Turbulent flows[END_REF]. In this region, the wallshear stress influences the flow. In the case of a zero-pressure gradient boundary layer over a flat plate, the advection terms are negligibly and a balance can be observed between the viscous and the turbulent shear stress. The flow inside the turbulent boundary layer can be rescaled based on the friction velocity but measuring directly the friction velocity in an experiment can be relatively difficult. To this end, the Clauser chart method [START_REF] Clauser | Turbulent boundary layers in adverse pressure gradients[END_REF] is used to estimate . By sketching + versus + (where + = ∕ ) on a semi-logarithmic scale along the -axis, the overlap region of the boundary layer is denoted as a straight line. This region is subsequent to the log law given by

( ) = 1 log + , (6)
where = 0.41 is the von Kármán constant = 5 is the constant used in [START_REF] Clauser | Turbulent boundary layers in adverse pressure gradients[END_REF]. The presence of the logarithmic region in an adverse pressure gradient boundary layer is still under discussion [START_REF] Monty | A parametric study of adverse pressure gradient turbulent boundary layers[END_REF], but it is known that a turbulent boundary under the effect of an adverse pressure-gradient evolves towards the zero-pressure-gradient state. [START_REF] Monty | A parametric study of adverse pressure gradient turbulent boundary layers[END_REF] pointed out that the skin friction, estimated from the Clauser chart method is matched with that determined from oil-interferometry for smooth adverse pressure gradient, at = 11860. In fact, for mild pressure gradients, that is when < 3, the skin friction obtained using both methods was satisfactory. However, for > 3, the difference between both approaches becomes significant, with up to 10% variation in . In equation ( 6), ( ) can be estimated directly from the experiments assuming that , , and are constants. Thus, the only unknown term is . Besides that, figure 5b shows the normalized boundary layer by and + , and the different regions associated to the viscous sub-layer + < 5, the buffer layer 5 < + < 30, the turbulent region 30 < + < 300, and the outer region. All profiles show a fairly good agreement between the simulation (for both models) and the experiment. Note that DDES performs yet better than RANS as shown in figure 5b, both in terms of thickness and mean velocity. Figure 6a shows the evolution of , calculated by the simulation and compared to that estimated using the Clauser method. Consequently, 15 points were extracted from the simulation data, covering locations along ⋆ in the range [-3, 0] across the top of the body surface. Figure 6b reports the growth of , ⋆ and along the -direction, where the integral parameters ( ⋆ , ) also increase until they reach their maximum value at = -0.4 , remaining constant until the body's end.

Global aerodynamic quantities

Aerodynamics coefficients such as the drag coefficient are key factors to diagnose the performances of road vehicles but they can be particularly difficult to measure in full-scale applications. There is a clear need for key performances indicators that can be upscaled to full-scale models. In fact, is known to be linked to other flow parameters, such as the recirculation region length and the base pressure drag. For this purpose, we introduced a summarised comparison between experimental and simulation results in this section and further carry on the analysis in the next subsection.

The drag coefficients is given by

= 2 2 ∞ , ( 7
)
where is the drag force and is the projected model frontal area ( = 0.112 m 2 ). The pressure coefficient writes

= 2 -∞ 2 ∞ , ( 8
)
where ∞ is the free-stream static pressure. The base drag was estimated by integrating the pressure coefficient as where ̄ is the time average estimated base-pressure coefficient while is the area covered by the pressure taps. Mean base-pressure distributions ̄ ( ⋆ = 0, ⋆ , ⋆ ) are shown in figure 7. Experimental results are comparable to the results of [START_REF] Volpe | Experimental characterization of the unsteady natural wake of the full-scale square back ahmed body: flow bi-stability and spectral analysis[END_REF] and Grandemange et al. (2013c) but several differences can be noted with both numerical simulations. Pressure recovery is in better agreement with the DDES than RANS but both simulations overestimate the pressure at the centre of the base. Note that the pressure coefficient calculated near the edges, where the boundary layer ends, are in relatively good agreement, hinting to the role of the wake in these discrepancies. The mean pressure coefficient at the base in the experiment is in relatively good agreement with the DDES (6% difference) but begins to deviate for the -SST (9% difference) (see table 2).

= 1 ∮ ̄ d , ( 9
)
The mean recirculation region is shown in figure 8, measured at the transverse mid plane and for each of the three methods. While experiments compare favourably with the mean recirculation region obtained from the DDES, the RANS approach is a far fetch. The recirculation length is defined by the mean separation line (defined by the isovalue ( , ) = 0), such that = max( | ⋆ =0 ), at the tip of the recirculation region. In table 2, Δ = ( . -. )∕ .

is computed as the difference between the simulation and the experiment. From these results, it is clear that the recirculation length predicted with the RANS approach is 30% larger than the experiment. Instead, only 7% separate the results between the DDES and the experiment. Note that the isosurface defined by ( ) = 0 will be used later in order to understand where the discrepancies between the model and the experiment come from. The differences between RANS and DDES carry similarly for the base pressure coefficient and the drag coefficient. RANS under estimates the drag coefficient by 7% and the base-pressure coefficient by 9% compared to the experiment. At the contrary, DDES is only 2% away from experimental measurements while the base pressure reaches 6% error. Clearly, a such discrepancy cannot come from the boundary layer itself since it was found to be nearly equivalent between RANS and DDES. The streamlines of the mean flow are shown in figure 8(a,c,e) where the recirculation region differs essentially at the bottom shear layer. In the RANS case, the flow near the end of the recirculation region shows a surprisingly large amount of negative vertical velocity. This in turn produces a weaker horizontal pressure gradient as shown in figure 8(f) for the RANS case, compared to the experiment and the DDES case (cf. figures 8(b,d)). An analysis will be carried out in the region to understand how the bottom shear layer interacts with the ground and the upper shear layer.

Wake topology

The Recirculation Region Interface (RRI), defined as the zero-horizontal velocity surface, such that is the three-dimensional mean separation line shown using a continuous black line in the two-dimensional plane in figure 8, and bounds the mass exchanges between the recirculation region and the outer flow. The RRI is a useful tool as it provides both a systematic geometric tool to measure, for instance, the recirculation length, or integrate the momentum equation to derive an equation relating the aerodynamic drag to the flow quantities in a budget. The streamlines depicted in figure 8 show that the RRI passes through the centre of the vortex of the mean wake and therefore accounts for the mean mass entrained within the wake [START_REF] Stella | Mass entrainment-based model for separating flows[END_REF]. The vortex-centre top location for the experiment and DDES model, along the upper shear layer is located at ⋆ =0.9 and ⋆ =0.82, while it is located at ⋆ =0.9 and ⋆ =0.15 for the lower shear layer. In contrast, the RANS simulation gives a distinct location, at ⋆ =0 for the lower shear layer and ⋆ =0.87 for the upper one both at ⋆ =1. This inconsistency in vortex location and dimension has a direct effect on the shape of the recirculation region, which can be clearly observed in figure 8 and thereby impacts the mass entrained together with drag.

( ) = ( ) | | | ( )=0 , ( 10
)

Shear-layer growth

As illustrated in figure 1, the RRI delimits the inner boundary of the shear layer and the respective development of both the upper and lower shear layers have to be related.

The separated shear-layer development is analysed based on the vorticity thickness , and momentum thickness . Following [START_REF] Dandois | Numerical simulation of active separation control by a synthetic jet[END_REF], and are computed as follows

= ( ) - ( ) ( ( , )∕ ) , ( 11
) and = ∫ ∞ 0 ( , ) - ( ) ( ( ) - ( ) 1 - ( , ) - ( ) ( ( ) - ( ) , ( 12
)
where ( ) and ( ) are the minimum, and maximum velocities respectively, along the streamwise direction. The shear layer is further divided into two parts (i.e. upper and lower) reconnecting at = . The shear-layer thicknesses linearly increase until a value of ∕ ≈ 0.8 for both, the experiment and DDES. For the RANS case, this growth stops near ∕ = 0.4. [START_REF] Bonnet | Collaborative testing of eddy structure identification methods in free turbulent shear flows[END_REF] suggest that the growth rate of the shear layer scales as = 0.17 According to equation ( 13), ∕ should be approximately equal to 0.2-0.23. In our case, the slope is approximately 0.34 near the body and 0.25 over the full extent of the near wake, which is in agreement with values expected for a free mixing layer, as reported for two-dimensional mixing layers [START_REF] Brown | On density effects and large structure in turbulent mixing layers[END_REF]) or the curved mixing layers [START_REF] Tenaud | Numerical simulation of the turbulent detached flow around a thick flat plate[END_REF].

Figure 9b shows the growth of the momentum thickness . The behaviour of is similar to but with a slope of 0.047. In their experiment of free mixing layers, [START_REF] Browand | The turbulent mixing layer: geometry of large vortices[END_REF] reported that the growth rate of the momentum thickness is well predicted by = 0.034 .

(

) 14
In reference to equation ( 14), ∕ should be equivalent to 0.041-0.046 which is consistent with the value of 0.047 expected for a free mixing layer. In addition, the rate of shear-layer growth ∕ is 4.6 for ∕ = 1 which is also consistent with the observations of a 2D mixing layer [START_REF] Bonnet | Collaborative testing of eddy structure identification methods in free turbulent shear flows[END_REF] where ∕ ≈ 4.8 was found. The growth of ∕ is 5.6 for ∕ = 1. With respect to their mixing layers, [START_REF] Castro | The structure of a turbulent shear layer bounding a separation region[END_REF] and [START_REF] Larchevêque | Large-eddy simulation of a compressible flow in a three-dimensional open cavity at high reynolds number[END_REF] observed that the ratio between the growth rates of the momentum and vorticity thicknesses has a value close to 5, which is comparable with that noticed in the present study at the lower part of the wake. Note that the symmetry between the upper and lower shear layer in the experiment can be explained by the rather large value of the ground clearance which has little effect for the present setup (see Grandemange et al. (2013c) for a detailed study on ground clearance effects).

The shear-layer growth reported for the experiment is in fairly good agreement with those observed from the DDES model, as shown in figure 9. At the contrary, RANS behaves differently. Indeed, both and are characterised by a lower shear-layer growth rate, where the ground proximity considerably impacts the shear-layer development. The separated shear layers act similarly to a free-shear layer, and their growth are dominated by the entrainment in the large region downstream Ahmed body model. However, this would have significant importance for the present study because the spreading of the free-shear layers is driven by the entrainment of external fluid [START_REF] Pope | Turbulent flows[END_REF]. In order to shed light on the sources of differences between DDES/experiments and RANS simulations, we turn to the analysis of the mass entrained into the wake to complement the results obtained from the shear-layer analysis.

Mass entrainment

Since the shear-layer growth drives the aerodynamic properties in the near wake, the mass transfer across the boundaries of the shear layer set the topology of the wake, and the RRI is used to diagnose the mass exchanges between the recirculation region and the outer flow. This topic is of fundamental importance for understanding and controlling the near wake because of two reasons. First, several studies have shown that one of the driving mechanisms of the near wake is the dynamic balance between the back flow and the amount of fluid that reenters the wake through the shear layer in the region of separation (Stella et al.,

| ̇ ⋆ | | ̇ ⋆ | 75% 75% 48%
Table 3 Mass flux integrated along the RRI interface (normalized by ∞ ).

2020

). Second, the results presented in §5.4 have shown that the size of the recirculation region is directly related to the growth rate of the shear layer. [START_REF] Gerrard | The mechanics of the formation region of vortices behind bluff bodies[END_REF] proposed a mechanism, which is the result of the equilibrium of a mass flux entering and leaving the recirculation region which sets the recirculation length. For the present case, the growth rate of the free-shear layers is associated with velocity entrainment along the shear layer. Therefore, it seems necessary to tie a link between shear-layer growth and the mass balance equation through the RRI which precisely follows the shear layer. The same principle of mass entrainment and recovery, between the shear and the recirculation regions, is expressed by [START_REF] Mair | Drag-reducing techniques for axi-symmetric bluff bodies[END_REF] for two-dimensional bluff bodies as

̇ ⋆ = -∫ ⋆ ( ⋆ ) cos ( ( ⋆ )) ⋆ = -∫ ⋆ ( ⋆ ) ⋆ (15)
where ̇ ⋆ is the mean mass flux, is the RRI boundary length, and is the angle between the normal of the surface and the -axis; while, ⋆ is the curvilinear abscissa at the RRI. In what follows, we only study a two-dimensional RRI obtained in the midplane at ⋆ = 0 for consistency with the experiment.

The evolution of the mass flux in the longitudinal direction, along the RRI, is shown in figure 10 where the mass flux is approximately symmetric. Therefore, the mass transfer according to its positive/negative values, can be separately studied. Table 3 reflects the integration of the mass entrained along the first range (negative mass flux), denoted by ̇ ⋆ and ̇ ⋆ , for the upper and lower parts of the RRI respectively. While the integration of the mass flux over the second range (positive mass flux); that is termed ̇ ⋆ and ̇ ⋆ , for the upper and lower sides of the RRI interface. Thus, the first interesting point that is reported in Table 3, is the symmetric behaviour between the upper and lower parts of the RRI interface for the experiment and DDES model. In contrast, RANS results are different. The entrainment feeding the recirculation comes from the upper part of the interface and explains the asymmetry of the recirculation bubble. Shear-layer and the mass-entrainment analyses show that the comparison between experiments and DDES holds fairly well but the RANS turbulence model clearly under predicts the entrainment at the bottom shear layer. In the following section, this comparison will be completed by linking this effect to the total drag. As shown in table 2, more than 80% of the drag comes from the base drag. In other words, to describe how the near wake interacts and affects the recirculation region topology, it is necessary to understand what drives the base pressure and therefore drag changes. The RANS recirculation presents a higher aspect ratio ⋆ , distinct from the DDES and the experiment where an increase of entrainment decreases ⋆ . This scenario is consistent with Roshko's models on the wake equilibrium and base pressure [START_REF] Roshko | On the wake and drag of bluff bodies[END_REF](Roshko, , 1993a,b) ,b) as well as more recently reported in [START_REF] Barros | Bluff body drag manipulation using pulsed jets and coanda effect[END_REF]. Next, we analyse the momentum budget to link the shape of the near wake to the drag through the RRI.

Momentum budget

In an analogy to the simplified two-dimensional wake, the streamwise mean momentum balance, through the closing separating line (i.e. the DSI) of the recirculation region was reevaluated by [START_REF] Balachandar | Properties of the mean recirculation region in the wakes of two-dimensional bluff bodies[END_REF] who showed that the base pressure coefficient is linked to the pressure, normal, and shear stresses along the DSI interface. However, the main issue in this study is that the unavailable pressure downstream the body, in the experiments, is impacting this discussion. To circumvent this issue, we will use the pressure estimation method recently proposed by [START_REF] Shanmughan | Optimal pressure reconstruction based on planar particle image velocimetry and sparse sensor measurements[END_REF] and use the RRI for the analysis instead. We now proceed with the analysis of the turbulent Reynolds stresses which are shown in figure 11 for the experiment and both numerical simulations in the mid plane of the wake. In order to compute the momentum budget, the RRI line is divided it into two parts; the upper part and the lower part as shown in figure 1. Figures (12c and to reach a maximum ⟨ ′2 ⟩ ⋆ = 0.04, then decreases until ⟨ ′2 ⟩ ⋆ = 0.02 at ∕ = 1. In figure 12d, the decreasing ⟨ ′ ′ ⟩ ⋆ along the RRI is similar between experiment and DDES. The shear stresses decreases until ∕ = 0.8 to reach a minimum ⟨ ′ ′ ⟩ ⋆ = -0.02, then it increases until it reaches ⟨ ′ ′ ⟩ ⋆ = 0 at ∕ = 1. These results show that the normal and shear stresses play an important role in the 0.5 < ∕ < 1 region. Figures also show the different order of magnitude of the shear stress so that max(⟨ ′ ′ ⟩ ⋆ ) is half the maximum of (⟨ ′2 ⟩ ⋆ ). The difference between DDES/experiments and RANS simulations clearly arises from the turbulent stresses which are fully modelled using an eddy-viscosity assumption for the RANS case.

For two-dimensional bluff body flows, the streamwise momentum balance on the RRI line writes

= ∫ - + 2 ⟨ ′2 ⟩ ⋆ -2 ⟨ ′ ′ ⟩ ⋆ d ⋆ -∫ + + 2 ⟨ ′2 ⟩ ⋆ -2⟨ ′ ′ ⟩ ⋆ d ⋆ , ( 16
)
where represents the body-base pressure, is the pressure distribution along the RRI line in addition to normal and shear components. In what follows, each term and their contribution in the budget along the upper and lower part of the RRI is analysed and compared with the experiment through a cumulative integral analysis of equation ( 16) along the direction. Figure 13b shows that the pressure term is shown along with the upper (blue) and lower (red) part of the RRI line for the DDES model where two regions appear. The first one is from 0 < ∕ < 0.41, in this region, the pressure terms for both parts are symmetric and provide a negative contribution to the drag (an actual net thrust) (Red-Blue>0). The second region starts from 0.41 < ∕ < 1, here we notice that the symmetry between the upper and lower part is broken, and the pressure term starts to play a role in increasing the pressure base. Figure 13c, the pressure term is shown along with the upper and lower part as well, for the -SST model. Two regions can be again identified. The first one is from 0 < ∕ < 0.59 where in this region, the pressure term for the two parts is asymmetric and provide a negative contribution to the drag (Red-Blue>0). The second region starts from 0.59 < ∕ < 1, where the pressure term starts to play a role in increasing the pressure base. The physics between the SSTand DDES is completely different, and the -SST shows a region of thrust for the lower part of the RRI while a region of drag in the upper RRI. Experimental results are reported in figure 13a and compare favorably with the DDES. Note that for = , the pressure contribution for the experiment falls very close to that obtained in the DDES and demonstrates the robustness of the pressure reconstruction algorithm [START_REF] Shanmughan | Optimal pressure reconstruction based on planar particle image velocimetry and sparse sensor measurements[END_REF].

In figure 14(a,c), the cumulative integral along the streamwise direction of the normal stress term is shown along with the upper (blue) and lower (red) part of the RRI line for the experiment and the DDES model respectively. The first region here is from 0 < ∕ < 0.31, in this region, the normal stress term for the two parts are symmetric and increase the pressure base (Red-Blue<0). The second region starts from 0.31 < ∕ < 1, here we noticed that the symmetry between the upper and lower part is broken too, and the normal stress term starts to play a role in decreasing the pressure base. In figure 14(e), the normal stress term is shown along with the upper and lower part as well, for the SSTmodel. The first region is from 0 < ∕ < 0.5, in this region, the normal stress term for the two parts are asymmetric and play a helper role (Red-Blue>0). The second region starts from 0.5 < ∕ < 1, here we noticed that this term starts to play a role in decreasing the pressure base. The physics between the experiment and DDES is the same, while both with SSTare completely different.

In figure 14(b) and 14(d), the shear stress term is shown along with the upper and lower part of the RRI line for the experiment and DDES model. The first region here is from 0 < ∕ < 0.185 where the shear stress terms for the two parts are symmetric and decrease the base pressure (Red-Blue<0). The second region starts from 0.185 < ∕ < 1.

Here we notice that the symmetry between the upper and lower part is still present, while the shear-stress term starts to play a role in increasing the pressure base. In figure 14(f), the shear stress term is shown along with the upper and lower part, for themodel SST. The first region is from 0 < ∕ < 0.05, in this region, the shear stress terms on each side of the RRI are asymmetric and decrease the drag (Red-Blue>0). The second region starts from 0.05 < ∕ < 1 where shear stress term starts to play a role in increasing the pressure base. The physics between the experiment and DDES remain the same, while the -SST over predicts the shear stresses while under predicting the normal stresses. It comes clear that the eddy-viscosity assumption is unable to accurately predict the turbulent stresses, leading to a deficit in mass entrainment through the bottom boundary layer, leading to the asymmetry of the wake in the RANS simulation and the over extent of the recirculation region.

Conclusions

The present study analyses the time-averaged wake properties of an Ahmed body at high Reynolds numbers comparing mean results from experiments and simulations at = 7.7 × 10 5 . Two different turbulence models are analysed and compared with experimental results using mass and momentum budgets performed along the recirculation region interface for the vertical plane positioned at the centre of the wake. Boundary-layer characteristics, shear-layer growth, and the near-wake mass and momentum budgets allow for quantifying the role of each part of the flow in driving drag. A comparison of the drag coefficients and characteristic length scales of the wake between the experiment and the delayed detached eddy simulation confirm a good agreement.

The shear-layer growth appears to follow the same scaling laws than previously observed for free-shear layers, and this, for both numerical simulations and the experiment, thereby preventing any definitive conclusions regarding a quantitative comparison. The analysis is then reconsidered in the framework of a mass and momentum budget along the recirculation region interface which allows for accessing both the mass entrained through the recirculation region (respectively detrained outside). This simple, yet effective analysis highlights that for RANS simulations, the amount of mass entrained by the bottom shear layer is clearly under predicted by a factor 4 when compared with the experiment. Then, the momentum budget along the recirculation interface allows for linking turbulent stresses with the pressure inside the recirculation region and the base of the body. To do so, the pressure reconstruction technique recently proposed by [START_REF] Shanmughan | Optimal pressure reconstruction based on planar particle image velocimetry and sparse sensor measurements[END_REF] is used to compute the mean pressure in the wake from experimental data and shows an excellent agreement with the DDES in the mid plane although the flow is essentially symmetric yet threedimensional. The momentum budget highlights that turbulent stresses are systematically under predicted in the RANS simulation and asymmetric, which is symptomatic of the calibration of the eddy viscosity. Whence, under predicting both the normal and shear stresses leads to the mass entrainment deficit reported in the mass budget. This in turn alters the growth of the shear layers, and thus, the pressure at the base predicted by the RANS simulation, leading to an unrealistic downwash of the recirculation region.

On the other hand, the delayed detached eddy simulation captures the same physical behaviour than the experiment and the contributions to the momentum budget all compare favourably with the experiment. Both the geometry of the near wake and the contribution of the different terms in the base-drag budget follow nearly identical trends. Overall, the recirculation region interface can be sought as a generic mean to decompose and diagnose drag in bluff-body wakes and understand the momentum balance in complex three-dimensional flows.

The tools developed in this study naturally extend to three-dimensional wakes as well as more realistic geometries. This particular aspect is currently under investigation in our group in the case of an Ahmed body and a realistic car geometry, representing interesting test cases to understand, model, and control the physical phenomena at the origin of drag for road vehicles.

Figure 1 :

 1 Figure 1: Sketch of the flow and wake geometry of an ahmed body showing the dividing streamline interface (DSI) in black with dashed lines and the recirculation region interface (RRI) with a (red-blue) dashed-dotted line. The colour background shows the vorticity non-dimensionalised with the free-stream velocity ∞ and the height of the body (adapted from Passaggia et al. (2021)).

Figure 2 :

 2 Figure 2: Sketch of the experimental setup.

Figure 3 :

 3 Figure 3: (a) Computational domain characterised by = 13.2 , = 1.9 and = 1.7 , the yellow dot shows the location where boundary-layer measurements were taken in the experiments. (b) The mesh specialized by three added boxes, box 1 covers the whole model while box 2 covers the under body area and box 3 is further refined in the near wake, downstream the body. (c) Mean recirculation region interface defined ( ( ) = 0) coloured with the mean pressure computed from the DDES and (d) an instantaneous snapshot of the DDES coloured by means of the instantaneous horizontal velocity.

Figure 4 :

 4 Figure 4: Comparison of the time-averaged boundary-layer profile measured at = -1.5 from the front part of the body for the simulations (continuous line) and experiment (symbols).

Figure 5 :

 5 Figure 5: Mean velocity profile along the -axis normalised in (a) by the free stream velocity and the boundary-layer thickness . (b) Same profile normalised using the friction velocity and the friction height = ∕ . The DDES is shown in red, the RANS in blue, and the experiment by means of symbols.

Figure 6 :

 6 Figure 6: (a) Evolution of the friction velocity, along the roof of the body at = 0, black dots represent experimental data. Red represents the DDES while blue reports results from RANS. (b) Streamwise evolution of the boundary-layer characteristic thicknesses , ⋆ and .

Figure 7 :

 7 Figure 7: Pressure distribution ̄ ⋆ ( = 0, ⋆ , ⋆ ) for (a) experiment, (b) DDES, and (c) -SST measured at the base of the body.

Figure 8 :

 8 Figure 8: Time-averaged horizontal velocity and pressure coefficient contours along the --plane at = 0 for (a,b) experiments, (c,d) DDES, and (e,f) -SST.

Figure 9 :

 9 Figure 9: (a) Vorticity thickness ( ) normalised by the height of the body against the streamwise direction normalised by the recirculation length for the upper and lower part. (b) Momentum thickness ( ) normalized by the height of the body for the upper and lower part. The DDES is shown in red, the RANS in blue, and the experiment by means of symbols.

Figure 10 :

 10 Figure 10: Mass entrainment comparison along the RRI where the DDES is shown in red, the RANS in blue, and the experiment by means of symbols.

Figure 11 :

 11 Figure 11: Time-averaged turbulent horizontal normal stress and shear stress contours along the --plane at = 0 for (a,b) experiments, (c,d) DDES, and (e,f) -SST.

Figure 12 :

 12 Figure 12: (a) Recirculation region interface. (b) pressure projected along the RRI, (c) normal stress, ⟨ ′2 ⟩ ⋆ along the upper part of the RRI interface. (d) Shear stress, ⟨ ′ ′ ⟩ ⋆ along the upper part of the RRI interface. The DDES is shown in red, the RANS in blue, and the experiment by means of symbols.

Figure 13 :

 13 Figure 13: Cumulative integral of the pressure term along the upper and lower part of the RRI interface for (a) experiment, (b) DDES, (c) -SST. Blue symbols show the upper part while red symbols show cumulative integral on the lower part of the RRI.

Figure 14 :

 14 Figure 14: Normal stress term along the upper and lower part of RRI interface for (a) experiment, (c) DDES, (d) -SST (blue:upper, red:lower). Shear stress term along the upper and lower part of RRI interface for (b) experiment, (d) DDES, (f) -SST.

Table 1

 1 Boundary-layer characteristics from experimental and simulation data.

		∞ ( ∕ )	∕	⋆ ∕	∕
	Experiment	45	0.069 0.0097 0.0076 1.29 5420 -0.26
	DDES	44.4	0.055	0.01	0.0072 1.37 5307 -0.22
	-SST	43.9	0.059	0.01	0.0072 1.36 5256	-0.1

Table 2 Table

 2

of drag coefficient , base-pressure drag coefficient ( ̄ ) and normalised recirculation region length ⋆ = ∕ .

  2017; Lorite-Diez, Jimenez-Gonzalez, Pastur, Martinez-Bazan and Cadot,

		Exp.	DDES	-SST
	̇ ⋆	0.05	0.06	0.05
	̇ ⋆	0.03	0.04	0.008
	̇ ⋆	-0.03	-0.03	-0.05
	̇ ⋆	-0.03	-0.03	-0.04

⋆ ⋆ ⋆ = -⋆ ⋆ -⋆ ⟨ ′ ′ ⟩ ⋆ -1 ⋆ ⋆ + ⋆ ⋆ ,(3)

Acknowledgement:

This work was supported by the French National Research Agency (ANR) through the COWAVE Project (ANR-17-CE22-0008-01). This work was granted access to the HPC resources of CINES/IDRIS under the allocations 2019-2a0129 and 2020-2a0129 made by GENCI.


domain: phys chim phys.cond

text:
Myroslav Holovko

Vojko Vlachy

Dung Di Caprio

D Di Caprio

On the contact conditions for the density and charge profiles in the theory of electrical double layer: From planar to spherical and cylindrical geometry

In this paper, from the Born-Green-Yvon equations of the liquid-state theory, we formulate two equivalent approaches for the calculation of the charge and density profiles of ionic fluids near charged nonplanar surfaces. The contact theorems for the density and charge profiles are formulated by the direct integration of this system of equations. Both theorems have nonlocal character. It is shown that the contact value of the density profile for noncharged surfaces is characterized by the bulk pressure and the surface tension. The obtained results are applied to spherical and and cylindrical surfaces. It is shown that the contact theorem for the density profile coincides with the recent results obtained by W.

The contact theorems are one of the few exact results in the theory of the electrical double layer. They establish the exact expressions for the contact values of the density and of the charge profiles in the interface between an electrolyte and a charged electrode. The first contact theorem was formulated more than fourty years ago by D. Henderson, L. Blum and J.L. Lebowitz [START_REF] Henderson | Some exact results and the application of the mean spherical approximation to charged hard spheres near a charged hard wall[END_REF][START_REF] Henderson | An exact formula for the contact value of the density profile of a system of charged hard spheres near a charged wall[END_REF] for the density profile in a planar electrical double layer formed by a primitive model of electrolyte and a uniformly charged planar hard electrode. According to their results, the contact value of the density profile is given by the sum of the bulk osmotic pressure P of the electrolyte and the Maxwell stress tensor contribution

ρ ct = α ρ α d α 2 = βP + β εE 2 8π ( 1
)
where ρ α (z) are the density distribution functions for ions of type α at a normal distance z from the wall, d α are the ion diameters, β = 1/(k B T ) with k B the Boltzmann constant, T the absolute temperature, ε is the dielectric constant of the solvent, where E is the electric field and εE 4π = q s is the surface charge density per unit area on the wall.

We note that in order to avoid the problem of electrostatic image effects between the electrode and the electrolyte, identical dielectric constant ε in the electrode and in the electrolyte was assumed. A little later, it was shown in [START_REF] Carnie | The statistical mechanics of the electrical double layer: Stress tensor and contact conditions[END_REF] that a formal derivation of the contact theorem for the density profile can be obtained by the direct integration of the Born-Yvon-Green (BGY) [START_REF] Yvon | La théorie statistique des fluides et l'équation d'état[END_REF][START_REF] Born | A General Kinetic Theory of Liquids[END_REF] equation between singlet and pair distribution functions for electrolytes in the presence of a charged electrode. In our previous papers [START_REF] Holovko | Contact conditions for the charge in the theory of the electrical double layer[END_REF][START_REF] Holovko | New results from the contact theorem for the charge profile for symmetric electrolytes[END_REF][START_REF] Holovko | On the contact condition for the charge profile in the theory of the electrical double layer for non symmetrical electrolytes[END_REF], in the framework, of a similar approach based on the integration of the BGY equation, the contact theorem for the charge profile was formulated for the planar electrical double layer.

In contrast to the contact theorem for the density profile, the contact value of the charge profile has a non local character and can be presented in the following form

q ct = α e α ρ α d α 2 = β α e 2 α ∞ dα/2 ρ α (z) ∂ψ(z) ∂z dz + β α e α P α (2)
where e α is the charge of the ion of type α, ψ(z) is the electrical potential defined as

ψ(z) = - 4π ε ∞ z q(z 1 )(z 1 -z)dz 1 (3)
where

q(z) = α e α ρ α (z) (4)
is the charge profile, P α is the bulk partial pressure for the ions of type α.

For the symmetrical electrolyte with d + = d -= d and e + = -e -= e, the first term in Eq. ( 2) equals zero due to the electroneutrality condition. As a result, the contact theorem Eq. ( 2) reduces to the more simple form

q d 2 = βe 2 ∞ d/2 ρ(z) ∂ψ(z) ∂z dz (5)
where ρ(z) is the density profile.

From the analysis of the effects of the ionic size on the properties of the diffuse double layer, using experimental and computer simulations results for symmetrical electrolytes for small charge at the wall, the semi empirical local expression for the contact value q(d/2), the so called Fawcett-Henderson-Boda conjecture [START_REF] Fawcett | A simple model for the diffuse double layer based on a generalized mean spherical approximation[END_REF][START_REF] Henderson | On a conjecture of Fawcett[END_REF] has been proposed. According to this expression

q d 2 = βe 2 1 χ P E ( 6
)
where χ is the inverse Debye length.

The connection between expressions ( 5) and ( 6) was the subject of discussions in [START_REF] Bhuiyan | Comment on "Contact conditions for the charge in the theory of the electrical double layer[END_REF][START_REF] Holovko | Response to "Comment on 'Contact conditions for the charge in the theory of the electrical double layer[END_REF].

In particular, it was shown that for a symmetrical electrolyte at small surface charge both expressions give similar results which are in good agreement with computer simulation data [START_REF] Bhuiyan | Comment on "Contact conditions for the charge in the theory of the electrical double layer[END_REF]. In [START_REF] Holovko | Response to "Comment on 'Contact conditions for the charge in the theory of the electrical double layer[END_REF], it has been established that expression (6) can be obtained from Eq. ( 5) as a result of two approximations. The first approximation is connected with the replacement in the integral in Eq. ( 5) of the density profile ρ(z) by its contact value Eq. (1). As a result

q ct d 2 = βe 2 ρ ct d 2 ψ d 2 (7)
where ψ(d/2) is the potential at the surface. The second approximation is connected with the use of the linear Gouy-Chapman theory [START_REF] Gouy | Sur la constitution de la charge électrique à la surface d'un électrolyte[END_REF][START_REF] Chapman | A contribution to the theory of electrocapillarity[END_REF] for the description of ψ(d/2).

The generalization of the contact theorems for charged non planar surfaces is one of the challenges in the theory of the electrical double layer. In particular, after the formulation of the contact theorem for the density profile in the case of the planar double layer, relation (1), significant efforts have been focused on the generalization of this contact theorem for the double layers with spherical, cylindrical and more complex non planar geometries. We note that according to the Henderson-Abraham-Barker approach [START_REF] Henderson | The Ornstein-Zernike equation for a fluid in contact with a surface[END_REF] the surface in the fluid can be considered as an additional particle in the infinite radius and infinite dilution limit. As a consequence, the density profile can be obtained as the pair distribution function between this particle representing the wall in the limit and a fluid particle. It was also shown in [START_REF] Henderson | Statistical mechanics of fluids at spherical structureless walls[END_REF][START_REF] Henderson | Statistical mechanics of fluid interfaces in cylindrical symmetry[END_REF] that in order to take into account the surface curvature effect, a surface term should be added to the pressure in expression [START_REF] Henderson | Some exact results and the application of the mean spherical approximation to charged hard spheres near a charged hard wall[END_REF]. This means that the excess chemical potential for the particle created can be presented in the form

µ ex i = P V i + γ i d 2 + R i S i (8)
where γ i (r) is the surface tension defined as the surface part of the chemical potential, V i is the volume accessible to the system, S i is the surface area of the curved wall. The subscript i = s or c denotes the spherical or cylindrical wall with the radius

R i = R s or R i = R c respectively.
We note that an expression similar to (8) was also used in the scaled particle theory (SPT) [START_REF] Reiss | Aspects of the statistical thermodynamics of real fluids[END_REF][START_REF] Reiss | Further development of scaled particle theory for rigid spheres: Application of the statistical thermodynamics of curved surfaces[END_REF][START_REF] Holovko | Fluids in random porous media: Scaled particle theory[END_REF] for the description of the chemical potential of the scaled particle in the hard sphere fluid. As for the SPT,

γ i (d/2 + R i ) can be expanded around (d/2)/(d/2 + R i ).
In a similar way, in the case of the hard sphere fluid near a spherical surface, Bryk et al. [START_REF] Bryk | Hard-sphere fluids in contact with curved substrates[END_REF] obtained the following expression for the contact value of the density profile near a non charged wall (E = 0)

ρ E=0 s d 2 + R s = ρ HS s d 2 + R s = βP HS - 9η 2 4π(1 -η) 3 (d/2) 3 (1 + η) d/2 R s + d/2 -η (d/2) 2 (R s + d/2) 2 (9)
where is ρ s is the total density distribution in the case of the spherical surface, where η = 1 6 πρd 3 is the packing fraction of hard sphere fluid. As usual, in the SPT approach [START_REF] Reiss | Aspects of the statistical thermodynamics of real fluids[END_REF][START_REF] Reiss | Further development of scaled particle theory for rigid spheres: Application of the statistical thermodynamics of curved surfaces[END_REF][START_REF] Holovko | Fluids in random porous media: Scaled particle theory[END_REF], the hard sphere pressure is given by the Percus-Yevick compressibility.

Comparably to [START_REF] Bryk | Hard-sphere fluids in contact with curved substrates[END_REF], the expression for the contact value of the density profile for the hard sphere fluid near a non charged cylindrical wall was obtained in the following form

ρ E=0 c d 2 + R c = ρ HS c d 2 + R c = βP HS - 9η 2 4π(1 -η) 3 (d/2) 3 1 2 (1 + η) d/2 R c + d/2 (10)
where is ρ c is the total density distribution in the case of the cylindrical surface,

We should note that in the case of ionic fluids near a non charged wall, the pressure P and the surface tension γ i (r) will include also the corresponding Coulombic contributions.

The second aspect of the generalization of the contact theorem eq. ( 8) for the non planar case is associated to detailing the second term in the expression [START_REF] Holovko | On the contact condition for the charge profile in the theory of the electrical double layer for non symmetrical electrolytes[END_REF]. So far, the activity in such direction was focused on the generalization of the Maxwell tensor for the double layer with charged spherical and cylindrical surfaces, charged membranes and charged surfaces with more complex geometries [START_REF] Vlachy | An alternative approach to the osmotic coefficient of polyelectrolyte solutions[END_REF][START_REF] Vlachy | On the virial equation for the osmotic pressure of linear polyelectrolytes[END_REF][START_REF] Wennerström | The cell model for polyelectrolyte systems. exact statistical mechanical relations, Monte Carlo simulations, and the Poisson-Boltzmann approximation[END_REF][START_REF] Trizac | Wigner-Seitz model of charged lamellar colloidal dispersions[END_REF][START_REF] Mallarino | The contact theorem for charged fluids: from planar to curved geometries[END_REF]. However, in these publications, the ionic sizes were not taken into account neither explicitly nor implicitly in the calculation of the Maxwell stress tensor. It was shown that for electrolytes with point ions, the Maxwell stress tensor contribution to the contact values of the density profile for walls with different geometries has a similar form as in the case of the planar wall.

The exact formulation of the contact theorem for the density profile for electric double layers in spherical and cylindrical geometries was obtained by Silvestre-Alcantara, Henderson

and Bhuiyan [START_REF] Silvestre-Alcantara | Contact condition for the density profiles in spherical and cylindrical double layers[END_REF]. According to their results, the contact theorem for the density profile for the surface of nonplanar geometry can be presented in the form

ρ i R i + d 2 = ρ E=0 i R i + d 2 -βE ∞ R i +d/2 q(r)dr (11)
where q(r) = α e α ρ α (r) is the charge profile,

ρ E=0 i (R i + d/2
) is the contact value of the density profile near the non charged surface. The hard sphere contribution to the contact value of the density profile is given by expressions ( 9) and [START_REF] Henderson | On a conjecture of Fawcett[END_REF] for spherical and cylindrical surfaces respectively. As we can see the second term in [START_REF] Bhuiyan | Comment on "Contact conditions for the charge in the theory of the electrical double layer[END_REF] has a nonlocal character and is more universal than the first term.

In [START_REF] Silvestre-Alcantara | Contact condition for the density profiles in spherical and cylindrical double layers[END_REF], it was shown that in the limit R i → ∞, the contact theorem (11) reduces to

ρ d 2 = βP + βE 2 ε 4π (12)
which differs from the expression (1) by the second term which has a missing factor 1/2 compared to expression [START_REF] Henderson | Some exact results and the application of the mean spherical approximation to charged hard spheres near a charged hard wall[END_REF]. The authors explained this difference by the fact that in the planar case one must assume a second parallel electrode of opposite charge at infinite separation from the first electrode.

In this note, we will show that the contact theorem for the density profile in the form [START_REF] Bhuiyan | Comment on "Contact conditions for the charge in the theory of the electrical double layer[END_REF] can be obtained from the BGY equation. Starting from the BGY equation, we will also obtain the exact relation for the charge profile. For convenience, we present the BGY relation in the following form

∇ 1 ρ α (1) = -β∇ 1 V α (1)ρ α (1) -βρ α (1) γ ρ γ (2)g αγ (12)∇ 1 u αγ (r 12 )d2 (13)
where 1 and 2 stand for the position of the corresponding particles, V α (1) = e α V coul α (1) +

V sh α (1) is the external potential, u αγ (r 12 ) = e α e γ u coul (r 12 ) + u sh αγ (r 12 ) is the interaction pair potential. The subscripts α and γ stand for ions having charges e α and e γ , the superscripts coul and sh indicate the electrostatic and non electrostatic short-range contributions. ρ α [START_REF] Henderson | Some exact results and the application of the mean spherical approximation to charged hard spheres near a charged hard wall[END_REF] is the density distribution of the particle of type α at point 1, g αγ [START_REF] Holovko | Response to "Comment on 'Contact conditions for the charge in the theory of the electrical double layer[END_REF] is the pair distribution function between particles of type α and γ.

We consider here two equivalent approaches, where in both cases, we start from the BGY equation. In the first approach, we start from the BGY equation in the form

∇ 1 ρ α (1) = -βe α ρ α (1)∇ 1 V coul (1) -βρ α (1)∇ 1 V sh α (1) -βρ α (1)∇ 1 W α (1) (14)
where

∇ 1 W α (1) = γ d2ρ γ (2)g αγ (12) e α e γ ∇ 1 u coul (r 12 ) + ∇ 1 u sh αγ (r 12 ) (15)
corresponds to the average force on a particle of type α at point 1.

In the second approach, in a similar way as in [START_REF] Carnie | The statistical mechanics of the electrical double layer: Stress tensor and contact conditions[END_REF][START_REF] Holovko | Contact conditions for the charge in the theory of the electrical double layer[END_REF], we introduce the standard mean electrostatic potential ψ(1) defined by

∇ 1 ψ(1) = ∇ 1 V coul (1) + γ e γ ρ γ (2)∇ 1 u coul (12)d2 (16)
In this case, as in [START_REF] Holovko | Contact conditions for the charge in the theory of the electrical double layer[END_REF], we can rewrite the BGY equation in the following form

∇ 1 ρ α (1) = -βe α ρ α (1)∇ 1 ψ(1) -βρ α (1)∇ 1 V sh α (1) -βρ α (1)∇ 1 W α (1) (17)
where

∇ 1 W α (1) = γ d2ρ γ (2) e α e γ h αγ (12)∇ 1 u coul (r 12 ) + g αγ (12)∇ 1 u sh αγ (r 12 ) ( 18
)
corresponds to the average force on a particle of type α at point 1 with the exception of the contribution from the mean electrical potential and h αγ (12) = g αγ (12) -1 is the pair correlation function. Now we will formulate the equations for the density and the charge profiles

ρ(1) = α ρ α (1) and q(1) = α e α ρ α (1). ( 19
)
To continue, we can use the BGY equation in the standard form ( 14) which includes the Coulomb interaction between wall and ions, V coul (1), or in the form [START_REF] Henderson | Statistical mechanics of fluid interfaces in cylindrical symmetry[END_REF] which includes the mean electrostatic potential ψ(1). We note that in both cases, we can omit the term related to the external potential V sh α (1) which we consider as equivalent to a hard wall. As a result from equation ( 14), we can formulate the following system of equations for the density and charge profiles

∇ 1 ρ α (1) = -β∇ 1 V coul (1)q(1) -β α ρ α (1)∇ 1 W α (1) (20)
∇ 1 q(1) = -β∇ 1 V coul (1)ρ(1) -β α e α ρ α (1)∇ 1 W α (1) (21)
where

ρ(1) = α e 2 α ρ α (1) (22)
After integration of relations ( 20) and ( 21), we obtain the contact theorem for the density and charge profiles generalizing the corresponding relations we obtained in [START_REF] Holovko | Contact conditions for the charge in the theory of the electrical double layer[END_REF] for the planar case. We note that the first term in relations ( 20) and ( 21), the Coulomb potential reads

V coul s (1) ∝ 1 εr Q s for spherical surfaces, and V coul c (1) ∝ - 1 ε Q c ln(r/r c ) for cylindrical surfaces,
where Q s and Q c are the charge of the corresponding surfaces. Since

∇ 1 V coul s (1) ∝ - 1 εr 2 Q s and ∇ 1 V coul c
(1) ∝ -1 εr Q c after integration of equations ( 20) and ( 21), we can present the generalization of the contact theorems for the density and charge profiles near a nonplanar wall. For simplicity, we assume d + = d -= d and we obtain the contact theorems in the following form

ρ i R i + d 2 = ρ E=0 i,α R i + d 2 -βE α ∞ R i +d/2 q i (r)dr (23)
q i R i + d 2 = α e α ρ E=0 i,α R i + d 2 -βE α ∞ R i +d/2 ρα (r)dr (24)
where ρ E=0 i,α (R i + d α /2) is the contact value of the density profile for the noncharged surface of type i (=c or s) and particle of type α, ρ E=0 i,α (R i + d α /2) is the partial contact value of the density profile for the ion of type α near non charged surface type i.

We note that the contact theorem [START_REF] Vlachy | On the virial equation for the osmotic pressure of linear polyelectrolytes[END_REF] for the density profile is identical to the relation [START_REF] Bhuiyan | Comment on "Contact conditions for the charge in the theory of the electrical double layer[END_REF] obtained by Silvester-Alcantara, Henderson and Bhuiyan [START_REF] Silvestre-Alcantara | Contact condition for the density profiles in spherical and cylindrical double layers[END_REF].

Starting from the BGY equation in the form [START_REF] Henderson | Statistical mechanics of fluid interfaces in cylindrical symmetry[END_REF], which includes the mean electrostatic potential ψ(1) instead of the Coulombic potential V coul (1), we can write the equations [START_REF] Holovko | Fluids in random porous media: Scaled particle theory[END_REF] and [START_REF] Bryk | Hard-sphere fluids in contact with curved substrates[END_REF] for the density and charge profiles in the following equivalent form

∇ 1 ρ α (1) = -βq(1)∇ 1 ψ(1) -β α ρ α (1)∇ 1 W α (1) (25)
∇ 1 q(1) = -β ρ(1)∇ 1 ψ(1) -β α e α ρ α (1)∇ 1 W α (1) (26)
Now after integration of these equations, we can present the contact theorems for the density and charge profiles near a nonplanar wall in in the following form

ρ i R i + d 2 = ρ E=0 i R i + d 2 -β ∞ R i +d/2 q i (r) ∂ψ(r) ∂r r δ i dr (27)
q i R i + d 2 = α e α ρ E=0 i,α R i + d 2 -β ∞ R i +d/2 ρ(r) ∂ψ(r) ∂r r δ i dr (28)
where δ i = 2 for spherical wall and δ i = 1 for cylindrical wall.

For charge symmetric electrolytes e + = -e -= e and the first term in the contact theorem for the charge profile in the form ( 24) or ( 28

and we have a simple nonlocal relation between the contact value of the charge profile and the density profile

q i R i + d 2 = -βe 2 E ∞ R i +d/2 ρ(r) dr (30)
or in the form

q i R i + d 2 = -βe 2 ∞ R i +d/2 ρ(r) ∂ψ ∂r r δ i dr (31)
In the limit R s → ∞ for the spherical surface or in the limit R c → ∞ for the cylindrical surface, the contact theorems for the charge and the density profiles reduce to the planar case. In particular in [START_REF] Silvestre-Alcantara | Contact condition for the density profiles in spherical and cylindrical double layers[END_REF], it was shown that due to the local electroneutrality condition

- ∞ d/2 q(z)dz = εE 4π , (32)
the last term in expression [START_REF] Vlachy | On the virial equation for the osmotic pressure of linear polyelectrolytes[END_REF] for the contact value of the density profile can be presented, in this limit, in the following form

-βE ∞ R i +d/2 q(r)dr ----→ R i →∞ -βE ∞ d/2 q(z)dz = β εE 2 4π (33)
which differs from the corresponding term in relation ( 1) by a factor 1/2.

In [START_REF] Silvestre-Alcantara | Contact condition for the density profiles in spherical and cylindrical double layers[END_REF], it was noted that the electroneutrality condition for the spherical and the cylindrical double layers has respectively the following forms

- ∞ Rs+d/2 q s (r)r 2 dr = R 2 s εE 2 4π , (34)
- ∞ Rc+d/2 q c (r)rdr = R c εE 2 4π . (35)
Due to this, for the nonplanar case, the contact theorem for the density profile cannot be presented in simple local forms as for the planar case.

We should note that the last term in the contact theorem in the form ( 27) due to the Poisson equation

∇ 2 1 ψ(1) = - 4π ε q(1) (36)
can be presented also in the form we see that the integral in (37) cannot give the simple Maxwell stress term as in the case of the planar wall [START_REF] Holovko | Contact conditions for the charge in the theory of the electrical double layer[END_REF].

-β ∞ R i +d/2 q i (r) ∂ψ ∂r r δ i dr = β ε 4π ∞ R i +d/2 ∇ 2 1 ψ i (1)∇ 1 ψ i (1)
Thus in this paper starting from the BGY equation between the singlet and pair distribution functions of ionic fluids near a noncharged planar surface, we formulated two equivalent approaches to calculate the charge and the density ionic profiles. The first one includes the direct Coulomb interaction between ions and surface. In the second approach, this direct ion-surface interaction is transformed through the Poisson equation into the mean electrical potential. The contact theorems for the density and charge profiles are then formulated for nonplanar surfaces by direct integration of this system of equations. The results are presented then applied to spherical and cylindrical surfaces. Both theorems have nonlocal character. It is shown that the contact value of the density profile for noncharged surfaces is characterized by bulk pressure and surface tension. The expressions of the contact theorem for the density profile are identical to the result recently formulated by Silvestre-Alcantara, Henderson and Bhuiyan [START_REF] Silvestre-Alcantara | Contact condition for the density profiles in spherical and cylindrical double layers[END_REF]. By using our previous result for the contact theorems for anisotropic fluids near a hard planar wall [START_REF] Holovko | Contact theorems for anisotropic fluids near a hard wall[END_REF], we can also generalize the obtained results

for the orientationally ordered charged dispersions case.

[29]

  ) vanishes. Moreover ρ(1) is proportional to the density profile ρ(1) = e 2 ρ(1)

Acknowledgments M. Holovko and D. di Caprio, acknowledge funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie grant agreement N˚734276. V. Vlachy acknowledges a support of the Slovenian Research Agency fund (ARRS) through project J1-1708.