Nathaniel T Baker
email: nathaniel.baker@onera.fr

Daniel Diaz

Didier Bailly

Laurent David

Jean-Claude Monnier

Jean-Claude Monnier

B Leclaire

J.-C Monnier

Mechanical vs. phenomenological formulations to determine mean aerodynamic drag from stereo-PIV wake measurements

Keywords: Control volume approach, Aircraft wake survey, Stereo-PIV

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

In a world increasingly concerned with energy and cost savings, it has become critical for aerodynamicists to accurately measure and predict the effects of drag acting on air vehicles. This task is all the more difficult, as drag may result from drastically different physical phenomena, ranging from viscous skin friction to shock waves, via the mere presence of lift. The mandatory first step towards improving the efficiency of an airframe is thus to accurately describe and measure these various drag components. The direct measurement of total mean drag acting on a model using force balances has been tried and tested, and is nowadays customarily found on industrial and research experimental platforms. Though reliable, this method only gives a global assessment of drag, and makes it impossible to distinguish between its various sources.

This drawback may be bypassed by expressing the total aerodynamic force in terms of the surrounding flow's properties, by means of a global balance of momentum. The main difficulty associated with the aforementioned method then comes from accurately evaluating the pressure and velocity fields throughout a control volume, which can sometimes be quite large. Within this context, Particle Image Velocimetry (PIV) appears to be a measurement method particularly well suited for this task, as it can yield velocity fields throughout large domains. And indeed, Van Oudheusden et al. [START_REF] Van Oudheusden | Non-intrusive load characterization of an airfoil using piv[END_REF] [START_REF] Van Oudheusden | Evaluation of integral forces and pressure fields from planar velocimetry data for incompressible and compressible flows[END_REF] demonstrated that the control volume approach fed exclusively with planar PIV data was viable when considering the mean aerodynamic loads acting on a two-dimensional (2D) airfoil. Since then, the continuous improvements made to PIV whether in terms of accessing the out of plane velocity component, or increasing spatial and temporal resolutions have made it possible to infer ever more precisely the relationship between the forces acting on an object, and its surrounding flow. For instance, Ragni et al. [START_REF] Ragni | Surface pressure and aerodynamic loads determination of a transonic airfoil based on particle image velocimetry[END_REF] [START_REF] Ragni | Nonintrusive aerodynamic loads analysis of an aircraft propeller blade[END_REF], computed the time and phase averaged loads acting on straightly flying and rotating airfoils respectively. The lift was computed by integrating the pressure coefficient along the airfoil's contour, which had been determined beforehand by integrating the Navier-Stokes equations fed with stereo-PIV data. Drag, on the other hand, was deduced from PIV measurements of the momentum deficit in the wake. These PIV based loads were reported to be within 10% of those measured using conventional methods (pressure taps and pitot probes). Alternately, Unal et al. [START_REF] Unal | Force prediction by piv imaging: a momentum-based approach[END_REF], Kurtulus et al. [START_REF] Kurtulus | Unsteady aerodynamic forces estimation on a square cylinder by tr-piv[END_REF], David et al. [START_REF] David | On the non-intrusive evaluation of fluid forces with the momentum equation approach[END_REF] and Villegas & Diez [START_REF] Villegas | Evaluation of unsteady pressure fields and forces in rotating airfoils from time-resolved piv[END_REF] took advantage of time resolved planar PIV data to examine the relationship between the instantaneous loads acting on 2D profiles and the unsteady surrounding flow, and studied the drag in terms of convective, pressure and turbulent contributions. In particular, these studies enabled to associate periodic features of the loads and vortex shedding in the wake. More recently, Terra et al. [START_REF] Sciacchitano | Aerodynamic drag of a transiting sphere by large-scale tomographicpiv[END_REF], [START_REF] Terra | Aerodynamic drag determination of a full-scale cyclist mannequin from large-scale ptv measurements[END_REF] have estimated the drag acting on bluff bodies, based on the investigation of their wakes with large scale tomographic PIV. Notably, the former study investigates the drag acting on a sphere towed across quiescent air, and highlights the competition between the convective and pressure contributions to drag, as a function of the distance between the investigation plane and the model.

The control volume approach in the form discussed above still entangles the different drag components resulting from irreversible losses of energy (i.e. viscous, form and wave drag, which add up to profile drag), to those resulting from the presence of lift (i.e. induced drag). To illustrate this point, let us consider the wake vortices that typically appear at the edges of threedimensional wings, which result from a sudden pressure discrepancy at the crossing of the wing-tip. This sharp pressure gradient ultimately drives a strong crossflow around the wingtip, whose magnitude can reach, in the present study, up to 20% of the streamwise velocity component. In the process of setting up wake vortices around its body, the airframe looses energy in the form of induced drag. It is however impossible to attribute the overall induced drag production to either inertia or pressure alone, since both effects contribute to the phenomenon. Within this context, the work of Jone's [START_REF] Jones | Measurement of profile drag by the pitottraverse method[END_REF] can be seen as a first attempt to extract the profile drag resulting from viscous and form drag only, by considering the momentum deficit in the wake of two-dimensional (2D) airfoils flying through an incompressible fluid. This approach was further extended by Maskell [START_REF] Maskell | Progress towards a method for the measurement of the components of the drag of a wing of finite span[END_REF] and Cummings et al. [START_REF] Cummings | Analysis of the elements of drag in three-dimensional viscous and inviscid flows[END_REF], who presented a way to determine lift induced drag. They did so by accounting for the presence of a transverse flow in the downstream plane, and showed that the induced drag in fact simply boiled down to the integration of the transverse kinetic energy at that same location. Finally, Destarac et al. [START_REF] Destarac | Drag/thrust analysis of jet-propelled transonic transport aircraft; definition of physical drag components[END_REF] and Kusunose et al. [START_REF] Kusunose | Wave drag extraction from profile drag based on a wake-integral method[END_REF] extended the latter methods to the compressible flow regime, where profile drag may also be generated by the emergence of shock waves upon the model. In such regimes, it is helpful to link profile drag to the production of entropy. In particular, the latter authors put forward a wave drag extraction method, based on the observation that the profile drag resulting from viscous shear layers is rotational (hence can be determined from Crocco's theorem), whereas the profile drag resulting from shock waves occurs in irrotationnal regions of the wake (hence wave drag can be calculated from the Rankine-Hugoniot shock jump formula).

The goal of the present work is to compare the performances of both approaches in an experimental context, where measurements stem from PIV. For the sake of simplicity, the framework adopted here will focus on the incompressible regime. We will start by reviewing the mechanical and phenomenological drag breakdown methods in section 2. The experimental rig used for the present investigation will then be described in section 3, while the implementation and robustness of both methods will be assessed in section 4. Global drag results will eventually be presented and discussed in section 5.

Theoretical background

Mechanical drag breakdown

Let us consider an aircraft model immersed in an upstream uniform flow U ∞ e x of density ρ and dynamic viscosity µ. The cornerstone of the method consists in performing a global balance of momentum over the control volume V of contour ∂V and outwards pointing normal n, which is fully encompassing a model of wing area S ref (cf. Fig. 1). Considering a steady mean flow and time averaged quantities, the mean aerodynamic load F acting on the model is traditionally found to be

F = - ∂V ρ (U • n) U dS - ∂V ρ u ⊗ u • n dS - ∂V P n dS + µ ∂V ∇U + (∇U) T • n dS, (1)
where U = U e x + V e y + W e z and P respectively represent the mean velocity and mean static pressure fields. The term u ⊗ u refers to the Reynolds stress tensor, with u = u e x + v e y + w e z being the turbulent velocity fluctuation. The effects of viscous friction along ∂V are considerably smaller than any other term (about two orders of magnitude here), hence can be safely dropped. Assuming that the lateral boundaries are sufficiently far from the model, there is no flux of momentum across them, and the global balance simplifies to evaluating the difference of momentum between the up-and downstream planes of the control volume S ∞ and S 1 respectively. This assumption also implies that the turbulent stresses will be non negligible only in the downstream plane S 1 . Further invoking the conservation of mass, and the fact that integrating P ∞ over the closed contour ∂V is zero, the total drag coefficient

C D = 2 F • e x /ρ ∞ S ref U 2
∞ can be expressed in terms of a surface integral over the downstream plane S 1 of the control volume, and reads

C D = 2 S ref S1 ρ ρ ∞ U U ∞ 1 - U U ∞ + P ∞ -P ρ ∞ U 2 ∞ - ρ ρ ∞ u u U 2 ∞ dS.
(

) 2
In the above, subscript '∞' refers to the uniform freestream quantities. Eq. 2 gives a mechanical description of the production of drag, by associating it to the properties of the flow (namely convection, pressure and turbulent stresses), which concurrently take part in all the various drag sources. The mechanical drag breakdown C Mec D ensues, and may be written as

C Mec D = C Conv D + C Press D + C Turb D , (3)
where

C Conv D = 2 S ref S1 ρ ρ ∞ U U ∞ 1 - U U ∞ dS, (4)
C Press D = 2 S ref S1 P ∞ -P ρ ∞ U 2 ∞ dS, (5)
and

C Turb D = - 2 S ref S1 ρ ρ ∞ u u U 2 ∞ dS. ( 6
)
According to this formulation, a sufficiently large plane S 1 must be considered, in order to compute a precise estimate of the total mean drag acting upon the model.

Phenomenological drag breakdown

Since the seminal work of Betz [START_REF] Betz | A method for the direct determination of wingsection drag[END_REF], numerous phenomenological drag breakdown methods have been put forward, many of which are reviewed and compared in Méheut & Bailly [START_REF] Méheut | Drag-breakdown methods from wake measurements[END_REF]. We focus in this paper on the phenomenological drag breakdown put forward by the latter authors, which gives a complete drag breakdown, whose integration is reduced to the sole wake. The cornerstone of Méheut and Bailly's phenomenological drag breakdown is Eq. 2. Unlike them however, we keep the Reynolds stress term, which is fully accessible using stereo-PIV data. (Neglecting turbulent stresses was required when wake data was measured with five-hole and Pitot probes, which only gave access to mean velocity fields.) Eq. 2 is then specialized to the aerodynamic case by expressing the ratio of fluid densities using the perfect gas constitutive law, in which the stagnation pressure P i and temperature T i are introduced. The flow is supposed to be isenthalpic (i.e. there are no heat sources), such that the stagnation temperature remains constant and equal to its freestream value everywhere. Finally, the downstream plane is assumed to be located sufficiently far away from the model, such that the quantities measured downstream differ only slightly from those upstream. Within this small perturbation framework, one can define P i = P i∞ + δP i (P i∞ being the freestream stagnation pressure) and U = (U ∞ + δu) e x + δv e y + δw e z , where the perturbations δP i /P i∞ , δu/U ∞ , δv/U ∞ and δw/U ∞ are all small compared to one. Owing to an asymptotic expansion up to second order, Méheut & Bailly showed that Eq. 2 could be re-written in the form

C D = 1 S ref w - 2 γM 2 ∞ δP i P i∞ + M 2 ∞ -1 δu 2 U 2 ∞ dS + 1 S ref S1 δv 2 + δw 2 U 2 ∞ dS - 2 S ref S1 ρ ρ ∞ u u U 2 ∞ dS, (7)
where γ represents the heat capacity ratio, and Ma ∞ = U ∞ / √ γrT ∞ is the freestream Mach number (r = 287 J/kg/K being the specific gas constant for air). The equation above now has a strong physical meaning. Indeed, the latter authors showed that the first integral on the right hand side of Eq. 7 translated the irreversible losses of energy, and was thus attributable to profile drag. By contrast, the second integral on the right hand side of Eq. 7 is proportional to the transverse kinetic energy that has emerged downstream of the model. As such, it is representative of the induced drag resulting from trailing vortices. In addition, the works of Maskell [START_REF] Maskell | Progress towards a method for the measurement of the components of the drag of a wing of finite span[END_REF] and Cummings et al. [START_REF] Cummings | Analysis of the elements of drag in three-dimensional viscous and inviscid flows[END_REF] showed that the integration domain of this second integral could be reduced to the wake only, by rewriting the transverse kinetic energy in terms of the streamwise vorticity component ω = (∇ × U) • e x , and the 2D mean streamfunction

ψ = -(∂ 2 yy ω + ∂ 2 zz ω) according to S1 δv 2 + δw 2 U 2 ∞ dS = w ψ ω U 2 ∞ dS. ( 8
)
As one may notice, the longitudinal velocity gradients must be negligible in the wake plane for Eq. 8 to be valid, hence restricts the wake plane to be evaluated sufficiently far away from the model (typically no less than a chord away from the model). Finally, turbulent stresses naturally vanish outside of the wake since they originate from shear. It is therefore possible to reduce the integration domain of the third integral on the right hand side of Eq. 7 to the wake only, as well. In the end, the phenomenological total drag breakdown C Phen D reads

C Phen D = C Prof D + C Ind D + C Turb D , (9)
where

C Prof D = 1 S ref w - 2 γM 2 ∞ δP i P i∞ + M 2 ∞ -1 δu 2 U 2 ∞ dS, (10)
C Ind D = 1 S ref w ψ ω U 2 ∞ dS, (11)
and

C Turb D = - 2 S ref w ρ ρ ∞ u u U 2 ∞ dS. ( 12
)
As already stated, the integration domain of the phenomenological breakdown ( 9) is reduced to the actual wake w, while the control volume approach (3) requires an integration over the entire downstream plane S 1 (with the exception of C Turb D being in fact the same in both formulations). This particular feature of the phenomenological method makes it quite appealing experimentally speaking, as it drastically reduces the span of the domain to capture. There is however a tradeoff with reducing the integration domain to the wake. Indeed, doing so focuses all the information scattered throughout space in a much smaller region, and especially in the viscous cores of the wake vortices. As a result, though the phenomenological method allows for smaller measurements domains, it requires a fine spatial resolution to accurately capture sharper velocity gradients.

Pressure calculation method

Inspecting Eqs. 5 and 10 above shows that the knowledge of both the static and stagnation pressure in the wake plane is required in order to evaluate drag using the mechanical and phenomenological approaches. Computing the pressure field from PIV velocity data, under the constraint of an incompressible flow is a long standing issue, which has received considerable attention in the past years [START_REF] Van Oudheusden | Piv-based pressure measurement[END_REF], [START_REF] Rival | Load-estimation techniques for unsteady incompressible flows[END_REF]. Note that in the more true to life case of a compressible flow, one may still compute the pressure field from velocity data, by invoking isentropic and isenthalpic relations, as put forward by Van Oudheusden et al. [START_REF] Van Oudheusden | Evaluation of integral forces and pressure fields from planar velocimetry data for incompressible and compressible flows[END_REF]. Owing to the the low value of the free-stream Mach number M a ∞ = 0.10, this particular study is restricted to the incompressible framework.

The pressure calculation strategy adopted here is based on the work of Jeon et al. [START_REF] Jeon | Sequential least-square reconstruction of instantaneous pressure field around a body from TR-PIV[END_REF], and relies on integrating the pressure gradient of the Navier-Stokes equations

-∇P = ρ [U • ∇U + ∇ • u ⊗ u ] -µ∆U. ( 13
)
The different terms composing the pressure gradient are calculated from the PIV measurements by spatial second order schemes. The pressure field is then obtained by minimizing a functional built on the difference of the pressure gradient based on the PIV measurements and the estimated pressure gradient, which is equivalent to solving a Poisson equation. The field is generally divided in sub-domains, based on the amplitude of the pressure gradient. The sequential pressure reconstruction is initiated from the outer domain, where the highest measurement accuracy is expected, hence where a reliable pressure reference can be taken. The pressure field in the other sub domains is then computed by imposing Dirichlet boundary conditions stemming from the previously computed outer domain.

Once the mean static pressure is reconstructed, the pressure integration constant P BC is adjusted so that the average of the reconstructed pressure computed along the contour of the PIV domain coincides with the average of the static pressure computed from the isentropic relation

P BC = P ∞ 1 + M 2 ∞ γ -1 2 1 - U 2 U 2 ∞ γ/(γ-1) (14)
at the same location. The local stagnation pressure in the downstream plane eventually ensues from applying the local isentropic and isenthalpic relations to the previously determined local static pressure, following

P i = P     1 + M 2 ∞ γ -1 2 1 - U 2 U 2 ∞ 1 + M 2 ∞ γ -1 2     -γ/(γ-1) . ( 15
)
3 Experimental apparatus and instrumentation

Wind-tunnel setup

This study was conducted on a model of the SAC-CON geometry, a generic flying wing configuration that was introduced within the framework of aircraft control and stability analyses [12] [26]. This particular configuration was chosen here as it is known to induce complex vortical patterns above its upper surface at medium and high angles of attack [START_REF] Roosenboom | Stereoscopic particle image velocimetry flowfield investigation of an unmanned combat air vehicle[END_REF] Three different angles of attack α = 9 o , 13 o and 18 o were investigated. The wake was measured one chord away from the model's wingtips (cf. Fig. 2). The maximum blockage, which was encountered for α = 18 o , did not exceed 2%, hence no particular correction was applied to the data.

Measurement protocol

The stereo-PIV setup revolved around a dual cavity Nd:YAG Quantel Evergreen laser, which emitted at wavelength 532 nm with an output energy of 200 mJ per laser head, and two LaVision Imager sCMOS cameras with a 5.5 Mpixels resolution. The PIV images were processed using 32 × 32pixels 2 interrogation windows, with an overlap of 50%, which resulted in velocity fields with a spatial resolution of 3 mm. Time averages were performed over 5000 consecutive snapshots acquired at 5 Hz in order to guarantee reliable statistics. Due to its large span, the wake was split into four 480 × 300 mm 2 overlapping PIV frames, with an overlap of 50 mm between adjacent frames. The cameras were fixed on a motor driven table that enabled a precise positioning of the different frames with respect to each other in space. In order to minimize the impact of edge effects when combining all frames together, the outer most frames were positioned such that they would contain the entire tip vortex. Finally, the longitudinal gradients present in Eq. 13 were discretized by considering two additional planes, respectively located at x = c ref ± ∆x, where ∆x = 30 mm. Subsequently to stereo-PIV measurements, the left half of the model's wake was surveyed using a five-hole probe. In order to reduce the acquisition time, the probe was mounted on a traversing system which swept across the wake at the constant speed of 5 mm/s. Data points were recorded every 2 mm, such that the characteristic traveling time τ p of the probe between two measurements is τ p = 0.4 s. By comparison, the characteristic time of the flow τ u is reasonably associated to the turnover time of the tip vortex following τ u = / U ⊥ , where is the size of the tip vortex, and U ⊥ the magnitude of the cross flow. Based on the data given in Fig. 3, τ u ∼ 0.03 s, which is one order of magnitude smaller than τ p . As a result, the five-hole probe experiences many turnover times during its travel time from one measuring point to another. It is therefore reasonable to assume that it experiences a statistically steady flow despite its motion.

Alongside the stereo-PIV and 5-hole probe measurements, the lift, drag and side forces, as well as the rolling, pitching and yawing moments were measured using an in-house force balance. The balance's dynamic range for the drag component was ±1730 N, with an uncertainty lower than 1.7 N.

General characterization of the wake

Figures 3 and4 present the overall structure of the flow for α = 9 o and 18 o respectively. The mean flow is given in Fig. 3, in which the contour plot gives the mean streamwise velocity component, while the vectors illustrate the velocity field in the plane normal to the stream. Since the wake is slightly asymmetrical, all graphs have been evaluated on the left half of the model, for the sake of a fair comparison. The right plots of figures 3 and 4 are in fact mirror images of the actual flow.

For α = 9 o , the wake features a small vortex core, which forms at the very tip of the model around |y| = 0.45 m. This primary vortex is accompanied by a narrow vortex sheet, which develops as the flow smoothly passes around the body of the model. At this setting, the level of turbulent fluctuations is negligible, and reaches only a couple of percents in the very core of the tip vortex. By contrast, the wake for α = 18 o , is dominated by a large vortex in the form of a kidney bean. This particular shape is typical of several unsteady vortices interacting with each other. As a matter of fact, Schütte et al. [START_REF] Schütte | Flow physics analyses of a generic unmanned combat aerial vehicle configuration[END_REF] showed that for α = 18 o , the tip vortex merges with a secondary vortex originating from the apex of the model. The unsteadiness of this process is confirmed by the maps of turbulent fluctuations displayed in Fig. 4, which shows a high level of normal turbulent stresses in the tip vortex region, reaching up to 10% of the free-stream kinetic energy.

In all cases, a strong cross flow is induced in the vicinity of the tip vortices, whose intensity reaches, in average up to 10% for α = 9 o and 20% for α = 18 o .

Experimental methods

Accuracy of the reconstructed pressure fields

The crux of the matter relies here in reconstructing the static and stagnation pressure fields in the downstream plane, from PIV measurements. In order to assess how reliable these reconstructed fields are, Figs. 5 and 6 compare streamwise profiles of PIV-reconstructed static and stagnation pressure to profiles stemming directly from 5-hole probe wake measurements. In both figures, the profiles are plotted along the left side of the model, as 5-hole probe measurements were available only there. These two figures strikingly show that the reconstructed static and stagnation pressure fields compare quite differently to their 5-hole probe counterparts. Indeed, the reconstructed static pressure (Fig. 5) is systematically lower than the direct probe measurements, while the reconstructed stagnation pressure (Fig. 6) appears to be much closer to the probe measurements, except for α = 18 o . Interestingly, the reconstructed stagnation pressure at low angles of attack appears to be quite accurate despite the fact that its computation uses a rather crude estimate of the static pressure field by virtue of Eq.15. This observation therefore suggests that the underestimation of the reconstructed static pressure is in fact mitigated by the velocity field, making the reconstructed stagnation pressure a more robust quantity for the purpose of wake integral methods. This is confirmed in Fig. 7, which compares two different methods of computing the stagnation pressure, based on Eq. 15. The first method consists in deducing P i from the static pressure, as done previously, while the second method consists in computing P i by imposing a uniform pressure throughout the wake plane. Interestingly, this figure shows that the global features of the stagnation pressure field can indeed be captured with the velocity field alone. However, in order to capture the stagnation pressure losses in the core of the tip vortex, it is preferable to compute P i from a prior computed, albeit coarse, static pressure field.

The discrepancy between the static and stagnation pressure reconstructed from PIV data, and their di- The solid line is associated to the direct five-hole probe measurements; The dashed line refers to P i computed from Eq. 15 fed with the prior computed static pressure P ; The dash-dotted line shows the result of computing P i with Eq. 15 in which an isobar wake plane at P ∞ is imposed.

rectly measured counterparts is better quantified by defining C p and C pi as

C P = S1 (P ∞ -P ) dy dz, (16)
and

C P i = S1 (P i∞ -P i ) dy dz, ( 17
)
and subsequently the relative differences P and P i

P = C P 5HP -C P PIV C P 5HP , (18)
and

P i = C P i 5HP -C P i PIV C P i 5HP , (19)
where superscripts "5HP" and "PIV" refer to integrals evaluated with 5-hole probe and PIV measurements respectively. Table 1 reports these ratios for all angles of attack investigated.

Wake identification

As accurate as the stagnation pressure estimate may be (whether directly measured or reconstructed), it still yields poor global results once integrated throughout the downstream plane, as stagnation pressure losses do not exactly vanish outside of the wake due to noise and uncertainties. It thus becomes necessary to extract the wake from its surrounding to enforce that stagnation pressure losses are confined there. The wake is identified by investigating the stagnation pressure distribution function, which presents a sharp discontinuity. A threshold value on the maximum admissible stagnation pressure loss is then imposed based on the value at which this discontinuity occurs (cf. the dotted line in Fig. 8). Fig. 9 shows the outer envelop of all the points lying beneath the aforementioned threshold, and confirms that the discontinuity in the stagnation pressure's distribution function indeed coincides with the physical boundaries of the wake. From here on, the integration of the stagnation pressure is limited to this boundary. The procedure above is illustrated here for α = 13 o , but works equally well at other angles of attack.

Sensitivity to spatial resolution

Tables 2 and3 show the sensitivity of the mechanical and phenomenological methods to the spatial resolution of the PIV data, which is given by the size of the interrogation windows (WS). Unsurprisingly, the two drag contributions that appear to be the most sensitive to the spatial resolution are those involving the reconstructed pressure field, namely C Press D and C Prof D , with a variation of forty to fifty drag counts between the highest and lowest resolutions (one drag count represents a 10 -4 variation of the drag coefficient). Furthermore, decreasing the spatial resolution affects the convective and turbulent contributions by about 10 drag counts in each case, as a result of C Conv D 's and C Turb D 's integrands missing the finer flow structures. Interestingly, the induced drag term is the least sensitive, and shows a variation of only 6 drag counts between the finest and coarsest cases, despite this term representing the velocity gradients in the plane perpendicular to the stream.

Sensitivity to longitudinal discretization

Earlier work focused on reconstructing a pressure field from planar PIV highlighted the necessity to evaluate the derivatives in the out of plane direction (cf. for instance [START_REF] Van Oudheusden | Piv-based pressure measurement[END_REF]) in order to get the most accurate predictions. The purpose of the following test is to assess the impact of the stream-wise gradients on the global drag predictions. To do so, we consider the three following schemes for the approximation of ∂ x u in Eq. 13, where for any function u(x, y, z):

                 ∂u ∂x = 0 + O(1) ∂u ∂x = u(c ref + ∆x) -u(c ref ) ∆x + O(∆x) ∂u ∂x = u(c ref + ∆x) -u(c ref -∆x) 2∆x + O(∆x 2 ) . (20)
The first scheme is a shorthand for neglecting longitudinal velocity gradients altogether. In the following, an interrogation window of size 32 × 32 pixels 2 is used. Practically speaking, ∆x = 30 mm is taken to be the spacing between adjacent PIV planes, while x = c ref refers to PIV measurements obtained in the middle plane located one chord away from the model's wingtips (cf. Fig. 2). Table 4 shows that in the case at hand, neglecting longitudinal velocity gradients in the pressure reconstruction algorithm has a somewhat larger effect on C Press D , than C Prof D , with a variation of up to twelve drag counts for the former vs. four for the latter. This again seems to indicate that the stagnation pressure is a more robust quantity than the plain static pressure. In any case, the sensitivity of either C Press D or C Prof D to the discretization scheme appears to be marginal, when compared against their respective sensitivity to the spatial resolution of the data, which suggests that it is the latter that should be given priority when computing loads from wake data.

Results and discussion

To begin with, table 5 gives the global drag coefficients measured by the force balance, which will be used as a reference throughout the remainder of the text. The relatively high measurement dispersion of 8% for α = 9 o comes from the low drag value occurring at this angle of attack, with respect to the balance's dynamic range. Tables 6 and7 display the global drag predictions and their breakdown, stemming from the mechanical (Eq. 3) and phenomenological (Eq. 9) approaches, as well as the relative difference between the global predictions and the balance measurements. Before computing the drag coefficients, the wake of the model holder was masked to remove its contribution to the total drag. Failing to do so would have resulted in an artificially inflated drag coefficient, since the model was directly mounted onto the force balance, at the tip of the holder.

According to tables 6 and 7, the mechanical and phenomenological methods yield global drag predictions that lie within 10% of the balance measurements, except for the mechanical prediction at α = 13, which is off by almost 20%. Furthermore, the mechanical breakdown appears to systematically underestimate the balance measured drag coefficients. In light of Fig. 5, this behavior likely results from the inaccurate reconstructed pressure from PIV measurements, which tends to miss the amplitude of the depression inside the tip vortex core. By comparison, the drag coefficients computed using the phenomenological decomposition systematically overestimate those measured with the force balance.

Table 6 Total mean drag predictions and their breakdown, as given by the mechanical approach. Computations are performed using the finest PIV resolution of 3 mm (32×32 pixels 2 interrogation windows). Values in parentheses are relative differences with respect to the balance measurements. Table 7 sheds some light on the relative competition between the physical processes at play. In particular, the phenomenological breakdown displays a reversal in the competition between the induced and the profile drag, with the former dominating at low angles of attacks (α = 9 and 13 o ), and the latter at α = 18 o . By contrast such a reversal does not transpire in the mechanical breakdown, for which the convective term always dominates. This observation may be physically interpreted in the light of Schütte et al.'s study [START_REF] Schütte | Flow physics analyses of a generic unmanned combat aerial vehicle configuration[END_REF], which shows that α = 18 o coincides with the tip vortex 9) integrated along z. For both cases, the red line represents the cumulative integral of all terms: the reading at y = 0 gives the total drag produced on one half of the wing. 9) integrated along z. For both cases, the red line represents the cumulative integral of all terms: the reading at y = 0 gives the total drag produced on one half of the wing.

Table 7 Total mean drag predictions and their breakdown, as given by the phenomenological approach. Computations are performed using the finest PIV resolution of 3 mm (32 × 32 pixels 2 interrogation windows). Values in parentheses are relative differences with respect to the balance measurements. being at the onset of bursting. From this point onward, the tip vortex progressively detaches from the upper surface of the model, therefore decreasing lift, hence also decreasing induced drag. Simultaneously, the resulting flow separation promotes pressure drag, which leads to an increase in profile drag. Such observations can prove quite valuable when optimizing an airframe for a particular flight regime, since tackling induced or profile drag requires radically different strategies. For instance, optimizing profile drag hints at reducing skin friction or delaying the flow separation, while reducing induced drag is usually tackled by preventing the formation of trailing vortices with wingtip devices. Inferring such strategies based on the mechanical decomposition alone would not be possible, since the competition between its different contributions may not automatically translate into physically interpretable drag sources.

α C Prof D C Ind D C Turb D C Phen D 9 o 0
The previous argument is illustrated with a more local point of view in Figs. 10 and11, which present spanwise profiles of mechanical and phenomenological drag, obtained by respectively integrating Eqs. 3 and 9 along the z coordinate. Here, the cases α = 9 and 18 o are considered. As before, all graphs are evaluated along the left half of the wing, the right parts of figures 10 and 11 being mirror images. For α = 9 o , the phenomenological decomposition (Fig. 10, right) concentrates drag in a narrow region of the wake which coincides with the location of the tip vortex's core. By contrast, the mechanical method (Fig. 10, left) associates drag for some part at the location of the tip vortex, but mostly along the model's centerline, where a deficit in the longitudinal velocity component appears, as mean kinetic energy has been redistributed into the cross flow. For α = 18 o however, the region of the wake mainly associated to the momentum deficit is clearly associated to the tip vortex. One must remain cautious when interpreting the wake in terms of drag production, since the latter only gives a picture of what the flow encountered along the model's surface, where drag was actually produced via pressure and skin friction. It can nonetheless yield some insight, especially in the presence of dominating flow structures such as here.

Finally, Figs. 10 and 11 show the cumulative integral of all the different contributions to the mechanical and phenomenological breakdown, where the value read on the model's centerline corresponds to the total drag produced by one half of the model. It should be observed that though the mechanical and the phenomenological methods translate in radically different spanwise profiles, they eventually integrate to similar values (within their respective accuracy), thus confirming that the two methods described in this paper are in fact two different points of view of the same phenomenon.

To put the current work in perspective with the existing litterature, the global drag predictions obtained here (whether using the mechanical or phenomenological method) appear to be at par with those published in earlier studies dealing with the same issue. Indeed, as far as wake surveys using the phenomenological decomposition fed with 5-hole probe measurements is concerned, Brune [START_REF] Brune | Quantitative low-speed wake surveys[END_REF] reports a relative difference better than 1% with respect to balance measurements on the drag prediction of a 1.80 m long rectangular wing, while Crowder et al. [START_REF] Crowder | Airplane flowfield measurements[END_REF] reports a 13% difference between balance measurements and the wake drag prediction of civil aircraft model flying at M a = 0.86. In this case, the main source of error is attributed to the method not accounting for flow compressibility effects.

Regarding load evaluation using the mechanical approach fed with PIV data, Ragni et al. reports drag values that lie within 10% of those obtained using a pitot wake rake behind a NACA0015 airfoil, while the drag prediction acting on a 2D square cylinder of Kurtulus et al. underestimates the expected value by 20% to 30%. In addtion, De Kat & Bleischwitz [START_REF] De Kat | Towards instantaneous lift and drag from stereo-piv wake measurements[END_REF] measured drag coefficients from a wake survey behind a flat plate, and reported drag coefficients within 30% of direct force measurements. The latter authors computed the total mean drag based on Eq. ( 2) fed with a flow field reconstructed from PIV measurements of a sub-region of the wake. The flow in the entire downstream plane S 1 was then inferred by invoking its symmetry with respect to the centerline, and its solenoidality. They thus attributed the relatively large discrepancies to the fact that they did not have access to the entire flow field. In that respect, the phenomenological method would perhaps be more suited, as it naturally vanishes outside of the wake. Based on the present work however, another aspect that may have had an impact in their results might be their assumption of an exactly symmetrical wake. Though reasonable in a first approximation, Fig. 9, and the cumulative values of Figs 10 and 11 which are systematically lower than half the total drag coefficient gives evidence that this may not always be exactly true.

Conclusion

This paper benchmarks the phenomenological drag breakdown method against the traditional control volume approach based on a global balance of momentum, using stereo-PIV data acquired in the complex wake of a finite-size wing. The former method expresses drag in terms of profile and induced drag, while the latter expresses drag in terms of the surrounding flow's properties. Until now, the phenomenological decomposition was numerically evaluated with five-hole probe measurements.

In the particular case at hand, we have shown that the mechanical and phenomenological methods perform similarly, as they yield integrated drag coefficients almost always within 10% of the balance measurements. Furthermore, we have illustrated how the phenomenological method can be used to physically interpret the drag sources at play, which can ultimately be used to optimize an airframe's geometry for a particular flight regime. By contrast, the mechanical approach does not allow one to distinguish between the different physical origins of drag, as it entangles the convective and pressure effects together.

As far as the implementation of the methods is concerned, this study highlights by comparing the PIVreconstructed static and stagnation pressure to direct five-hole probe measurements, that the latter is a more robust quantity than the former. This behavior likely results from the fact that the stagnation pressure in fact mitigates information about the reconstructed static pressure with information about the PIV-measured velocity field. In that respect, the phenomenological method appears to be particularly well suited to the investigation of aircraft wakes containing intense trailing edge vortices, as it enables one to reduce the integration domain based on physical arguments alone.

Fig. 1

 1 Fig.1Control volume of interest. The region enclosed in the dashed loop is referred to as the wake w, which is a subset of the downstream plane S 1 . Outside of this region, the profile drag, induced drag and turbulent contributions have completely vanished, unlike those associated to the convective and pressure terms.

Fig. 2

 2 Fig. 2 Overview of the stereo-PIV setup: (i) Quantel Evergreen laser; (ii) Laser sheet. In order to discretize the axial gradients present in equation (13), three PIV planes have been subsequently shone at locations x = c ref (solid line), and x = c ref ± ∆x (dashed lines); (iii) 5.5 Mpixels LaVision Imager sCMOS cameras; (iv) Wake plane split into four 480 × 300 mm 2 frames, with a 50 mm overlap.

  [START_REF] Schütte | Flow physics analyses of a generic unmanned combat aerial vehicle configuration[END_REF]. It is therefore an excellent case study to try and test wake integral methods. The particular model used here has a wingspan b = 1 m, a reference chord c ref = 0.31 m and a reference area S ref = 0.3253 m 2 . It was placed inside the L1 wind-tunnel of Lille's ONERA center, whose dodecagonal test-section is 2.4 m wide. The freestream velocity was fixed at U ∞ = 35 m/s throughout the campaign, which corresponds to a free-stream Reynolds number Re ∞ = ρ c ref U ∞ /µ = 1.1 10 6 , and a free-stream Mach number Ma ∞ = U ∞ / √ γrT ∞ = 0.10.

Fig. 3 2 ∞Fig. 4

 324 Fig. 3 Mean flow field in the reference plane x = c ref for α = 9 o (left) and α = 18 o (right). The colormap outlines the longitudinal velocity component U = U•e x , while the vectors represent the transverse flow in the wake plane U ⊥ = V e y +W e z .The mean ratio of the transverse to longitudinal velocity U ⊥ /U computed throughout the wake is indicated in the boxes.

αFig. 5

 5 Fig.5Spanwise profiles of static pressure integrated vertically. Solid lines refer to direct pressure measurements with a 5-Hole probe. Dashed lines refer to the reconstructed static pressure based on PIV measurements fed into Eq. 13.

αFig. 6 αFig. 7

 67 Fig. 6 Spanwise profiles of stagnation pressure integrated vertically. Solid lines refer to direct pressure measurements with a 5-Hole probe. Dashed lines refer to the reconstructed stagnation pressure based on PIV measurements fed into Eqs. 13 and 15.

Fig. 8

 8 Fig.8Distribution function of P i /P i∞ for α = 13 o . The dotted line locates the discontinuity in the histogram, whose argument is taken as the threshold value between the actual wake and its surrounding.

Fig. 9

 9 Fig.9Wake contour for α = 13 o , as identified by the threshold value found in Fig.8. This figure showcases the importance of cropping the wake's surrounding before integration, as the latter introduces artificial stagnation pressure losses, which have no physical meaning.

Fig. 10

 10 Fig. 10 Spanwise profiles of drag contributions for α = 9 o . (Left) mechanical decomposition: each black line represents a particular term of Eq. (3) integrated along z; (Right) phenomenological decomposition: each black line represents a particular term of Eq. (9) integrated along z. For both cases, the red line represents the cumulative integral of all terms: the reading at y = 0 gives the total drag produced on one half of the wing.

Fig. 11

 11 Fig. 11 Spanwise profiles of drag contributions for α = 18 o . (Left) mechanical decomposition: each black line represents a particular term of Eq. (3) integrated along z; (Right) phenomenological decomposition: each black line represents a particular term of Eq. (9) integrated along z. For both cases, the red line represents the cumulative integral of all terms: the reading at y = 0 gives the total drag produced on one half of the wing.

Table 1

 1 Relative differences between C P and C P i computed with 5-hole probe data and PIV measurements.

	α	P	P i
	9 o	0.85	0.04
	13 o	0.45	0.07
	18 o	0.18	0.17

Table 2

 2 Sensitivity of the mechanical decomposition to spatial resolution, for the α = 13 o case. Computations are performed using PIV interrogation windows of increasing sizes (WS).

	WS [pixels 2 ]	C Conv D	C Press D	C Turb D	C Mec D
	32 × 32	0.0294	0.0164	-0.0016	0.0442
	64 × 64	0.0293	0.0157	-0.0012	0.0438
	128 × 128	0.0283	0.0128	-0.0007	0.0404

Table 3

 3 Sensitivity of the phenomenological decomposition to spatial resolution, for the α = 13 o case. Computations are performed using PIV interrogation windows of increasing sizes (WS).

	WS [pixels 2 ]	C Prof D	C Ind D	C Turb D	C Phen D
	32 × 32	0.0205	0.0390	-0.0016	0.0579
	64 × 64	0.0196	0.0391	-0.0012	0.0575
	128 × 128	0.0151	0.0384	-0.0007	0.0528

Table 4

 4 Sensitivity of C PressD and C Prof Dto the order of the discretization scheme used in Eq. (13), for the α = 13 o case. Results are obtained with PIV interrogation windows of size 32 × 32 pixels 2 .

	Scheme order	C Press D	C Prof D
	0	0.0152	0.0201
	1	0.0158	0.0207
	2	0.0164	0.0205

Table 5

 5 Reference drag coefficients measured using a force balance, for the angles of attack α = 9 o , 13 o and 18 o .

		C Bal D
	α = 9 o	α = 13 o	α = 18 o
	0.0233 ± 0.0018	0.0543 ± 0.0016	0.171 ± 0.0034
	( ± 8%)	( ± 3%)	( ± 2%)

Acknowledgements The authors acknowledge the financial support of the Agence National pour la Recherche and the DGA under the grant ANR-16-ASTR-0005-01, as well as the CPER-FEDER of the Hauts-de-France region.