Drag decomposition of a subsonic wing via a far-field, exergy-based method

This paper focuses on the aerodynamic analysis and drag decomposition of an unpowered, low aspect ratio wing, using a far-field, exergy-based method. As opposed to traditional drag accounting methods, exergy balance provides insights into the amount of energy that can be potentially recovered off the body’s wake, which further translates into potential efficiency gains of the integrated engine-wing system. In this study, a far-field exergy balance method was used to determine the total drag of a three-dimensional wing. The far-field drag prediction was verified against near-field calculations. In addition, drag decomposition using exergetic terms was conducted to identify drag components that contain possibly recoverable energy. Such analysis can be subsequently used to educate the integration of a propulsion system to exploit the potentially recoverable wake energy and deliver an integrated engine-wing system with enhanced installed efficiency. The present methodology is a major step ahead in the application of far-field methods on three-dimensional wake domains and can potentially become a major enabler for optimal propulsion integration for future, novel aircraft-engine configurations.


Introduction
The global climate challenge drives the industry towards new means of aircraft propulsion aimed at eliminating CO 2 emissions. 1 However, the massive energy intensity of sustaining flight at both transonic as well as lower Mach number regimes increase the difficulty of the challenge. Therefore, it is even more important than before to seek means to reduce this energy intensity, and recuperation of energy that was initially wasted becomes a duty. It is expected that part of the envisaged improvements will be enabled via a more closely coupled integration between the airframe and the propulsion system. 2 Aircraft's wake recuperation, using more efficient installation of the propulsion system, was previously identified as a possible approach to increase the aerodynamic efficiency of the integrated system as shown by Smith. 3 However, conventional methods currently used to characterise the efficiency of certain installed configurations offer very limited insights into the potential of the propulsion system to recuperate part of the wake energy. A typical example of wake recuperation is the boundary layer ingestion concept a variation of which was previously analysed by Samuelsson and Gronstedt. 4 These systems take advantage of low flow momentum regions of the fuselage and/or wing boundary layer that is ingested by the propulsors, to enhance their propulsive efficiency. 5 A further concept is associated with the potential energy recovery of liftinduced vortex and it is represented by propellers mounted at the wingtips. 6 Even though this concept has been repeatedly studied by numerous authors in the past 7-10 the mechanisms that control the aero-propulsive benefits are still poorly understood. Early studies during the 60s and the 80s, reported potential benefits of up to 25% due to induced drag reduction and/or propeller propulsive efficiency increase. [7][8][9] Mechanical and aircraft control issues however created barriers to a further investigation. More recently, novel aircraft designs such as hybrid-electric or distributed fully electric propulsion concepts, have restored the interest of the research community in wingtip mounted propeller configurations, as an enabler for closer and more efficient integration between the engine and the aircraft. 10 Aircraft drag and thrust accounting, has been a challenging subject in the past both computationally as well as experimentally. 11,12 More specifically, in closely coupled designs, where the system's efficiency is highly dependent on the integration between the engine and the aircraft, the complexity of the task increases as there is no clear distinction between thrust and drag terms. 13 For unpowered conditions, a drag decomposition into its components is required to determine the impact of induced flow phenomena on the aerodynamic body performance. 14 The aim of such drag decomposition is to educate the aerodynamic integration of the propulsion system onto the aircraft as previously shown by Goldberg et al. 15 Additionally, for power-on conditions, a clear distinction of the aeropropulsive benefits and penalties, due to installation effects, is still required. 16,17 Conventional near-field methods, for the aero-propulsive assessment of highly integrated systems, despite their widespread use, provide poor insights into such drag and thrust breakdown. 18,19 The development of a commonly accepted aerothermo-propulsive methodology that addresses this issue, has attracted much interest by the research community in the past decades. 13 An extensive literature review of various drag analysis methods was presented by van Dam et al. 13 A series of numerical studies, aiming to provide aircraft drag prediction and mitigation, were also well documented by Deconinck et al., 20 providing detailed analysis of the benefits and drawbacks of each available method. More recently, the formulations of far-field drag decomposition methods and their applications for aircraft designs under poweroff and on conditions, have been extensively reviewed by Fan et al. 21 Considering the requirements of tightly coupled configurations, the most promising and recent formulations, have been developed and reported by Drela 22 and Arntz et al. 23 Drela 22 proposed a formulation which relied on mechanical power and kinetic energy flow analysis against the traditional forces and momentum-based methods. This methodology was validated on certain two-dimensional airfoil test cases with embedded wake ingestion technologies. A clear identification and quantification of the power sources, sinks and their interactions that influence the flight power requirements was provided. Sanders et al., 24 showed an application of Drela's formulation on a 2D boundary layer ingesting configuration under transonic conditions, with and without propulsion system in place. Arntz et al., 23 developed a formulation, for total energy management accounting also the thermal effects, which was the key difference from Drela's method. Arntz's mathematical formulation uses a combination of the first and second laws of thermodynamics enabling the so-called, exergy-based analysis. It requires no distinction between thrust and drag, allowing an integrated analysis, suitable for aircraft architectures that benefit from wake recuperation systems. The total supplied propulsive exergy is determined by the exergy balance. This is represented by the net propulsive power absorbed by the system and the losses, which are separated into reversible and irreversible phenomena. For the cases where no propulsive force is present (unpowered configurations), the net forces represent the total drag of the aircraft. The decomposition of the aircraft's drag, into its physical components, is determined by the energetic status of the wake. The interaction of the aerodynamic body with the flow generates a status of non-equilibrium due to velocity, temperature and pressure perturbations within the systems wake. These flow gradients represent the part of the losses that can be theoretically converted into mechanical work. Far-field drag decomposition provides critical information, for the aerodynamic analysis, related to the sources and nature of the overall system's drag. Arntz et al. 25,26 validated their mathematical formulation on unpowered and powered configurations. Initially, the Common Research Model was tested at transonic and cruise conditions. 25 Under power-on conditions the investigation focused on the boundary layer ingestion concept via a two-dimensional test case under Mach numbers of 0.2, 0.5 and 0.7, respectively. 26,27 These previous studies concluded that the concept of exergy provides a potential aero-thermopropulsive performance assessment methodology pertinent to boundary layer ingestion systems. Aguirre et al., 28 showed a further application of exergy balance on a two-dimensional NACA0012 airfoil and a threedimensional, unpowered wing. In this work, the benefit of the far-field method to provide richer insights into the involved physics was shown. The exergetic terms and total far-field drag however were only calculated at the trailing-edge of the bodies and at a single wake point without any further downstream Trefftz plane wake variations which would have examined the consistency of the exergy balance on the total drag calculations within the wake.
For the purposes of a commonly accepted aerothermo-propulsive methodology, suitable for assessment of unpowered aerodynamic bodies and closely coupled aircraft-engine designs, the application of farfield based formulations is yet to be widely shown. An attractive example for the application of such method is a flow-field dominated by the presence of the wingtip vortex which was previously identified as a flow mechanism that carries a significant amount of potentially recoverable energy. Wingtip vortex analysis has previously attracted attention by both academia and industry due to their critical impact on the overall wing drag performance. [29][30][31][32][33][34][35][36][37][38][39][40][41] In this work, a far-field, exergybased drag decomposition method, is demonstrated on a three-dimensional, low-subsonic, unpowered, rectangular NACA0012, wing configuration. The flow-field predictions are validated against experimental data by Chow et al. 42,43 The objective of the present work is to offer a detailed view of the physical phenomena linked to the aerodynamic drag of the wing and also to provide information about the potentially recoverable energy contained within the wake. The ultimate goal of the work is to enable the use of far-field wake analysis methods to determine the optimum propulsion system installation in future, novel aircraft configurations.

Exergy balance method
Exergy balance method was previously reported by Arntz et al. 23 The mathematical formulation relies on the first and second law of thermodynamics, momentum relation and mass conservation. This is restricted to mean steady flows based on the Boussinesq's hypothesis whereby to RANS, linear eddy viscosity turbulence models.
A control volume enclosing the wing S A and the outer surfaces S 0 is defined. A Trefftz plane on which wake parameters are calculated is defined at a certain location downstream of the wing's trailing-edge ( Figure  1). The axial position of the Trefftz plane is variable to enable the representation of the wake properties variation as a function of the downstream axial distance.
Following the definition of exergy, as the thermodynamic property that represents the maximum part of the energy that is theoretically fully convertible into mechanical work, one reads for an open system: where ε is the specific flow exergy, h i specific total enthalpy, s entropy and T denotes the static temperature. The above equation, accounts for no gravitational potential energy and air has been assumed as a perfect gas. The time-average spatial rate of change of exergy can be then expressed as follows: where ρ is density and V is the local velocity vector.
Integrating equation (2) within the control volume, using the divergence theorem, momentum relation, first and second laws of thermodynamics, the general exergy balance formulation for aero-thermo-propulsive performance assessment of integrated power-on configurations is derived. The mechanical and thermal power supplies of the propulsion system are associated with the exergy losses, entropy generation and power absorbed by the system to propel the aircraft. For unpowered conditions, an alternative form of exergy balance is available, representing a power-drag decomposition formula. The total drag of the fixed body is associated with exergy outflows and sinks. The exergy balance is then expressed as where D indicates the total drag, U ∞ the free-stream velocity, _ E tot total exergy outflows and _ A tot the total anergy. Each term is decomposed into mechanical (_ ε m ) and thermal (_ ε th ) outflows, viscous dissipation ( _ A Φ Þ, thermal mixing ( _ A =Τ ) and wave anergy ( _ A w ) respectively as follows: The exergy terms are surface integrals while the anergy ones are volume integrals. Each term in equation (4) is further decomposed as follows: Rate of mechanical exergy outflow ð_ ε m Þ: ÞðV Á nÞds : transverse kinetic energy deposition rateassociated with the trailing vortex system and represents the part of the induced drag that theoretically can be recovered.
A. Rate of thermal exergy outflow (_ ε th Þ: ρδsðV Á nÞds: rate of thermal anergy outflows denotes entropy. The three terms account for the maximum thermal power that can be recovered from the wake. As shown by Arntz et al. 23 thermal exergy outflow represents the thermo-compressible exergy available into thermal exergy rate at constant volume and thermal exergy rate linked with volume change. The thermal exergy outflow is then expressed mathematically via the following (for more details of the derivation see the detailed work of Arntz et al. 23 ): B. Rate of an energy generation due to dissipation ( _ A Φ ): where Φ eff , is the effective dissipation rate denotes an irreversible process linked to the available mechanical work in the wake due to the non-kinetic equilibrium and acts towards a homogenous velocity and pressure field. In addition, the dissipation process results to an increase in local temperature and therefore a gain in thermal exergy. Thus, the rate of change of anergy due to dissipation can be further decomposed as follows: C. Rate of anergy generation due to thermal mixing ( _ A =T ): where k eff , is the effective thermal conductivity (k eff ¼ C p ðμ=P r þ μ t =P r t Þ and the turbulent Prandtl number was chosen as a constant value of 0.9. _ A =T , denotes an irreversible process that destroys the available thermal work that exist in the wake due to the non-thermal equilibrium and acts towards a homogenous temperature field. D. Rate of anergy generation due to shock waves ( _ A w ): whereS w , is the sub-volume in which the irreversible process of shock waves occurs. _ A w returns the entropy risen due to this phenomenon. Since a low-subsonic test case is examined, the term is deactivated and will not be further investigated.
The exergy-based decomposition method is verified against the conventional, near-field drag (D nf ) calculation which considers the pressure and viscous forces as follows: where A ref , is the reference area of the aerodynamic body under study. The exergy balance terms are also nondimensionalised to express power counts (pc), here called far-field drag CD ff :  42 was respected in the directions of y=C (hight of 0.66 chords) and z=C (width of 1 chord) axes (see Figure 2). The inlet surface boundary was moved three chords upstream from its originally reported position (0.15 chords from the wing's leading-edge) for stability purposes. The axial distance from the wing's trailing-edge to the outlet boundary was also extended to four chords (from its original distance of 0.66 chords) to avoid flow circulation at the outlet and to enable flow analysis along a long distance. Following that and for the purposes of exergy-based analysis, a domain size study was conducted (gradually increasing the size of y=C and z=C axes) to ensure all the flow gradients have been integrated and blockage effects have been removed.
Steady-state, RANS based simulations were performed using ANSYS/Fluent 19.2 with the pressure-based solver and ideal-gas assumption. The governing equations system was solved under the scheme of the coupled algorithm using the spatial discretisation of the least squares cellbased method. Second order accuracy was implemented for pressure and second order upwind for density, momentum, energy and turbulence equations. The underrelaxation pseudo-transient method was used to improve numerical stability. The wing surface was modelled as a non-slip wall, whereas for the wind-tunnel walls, the zero-shear stress option (slip condition) was used. For the inlet and outlet planes, pressure-inlet and pressure-outlet boundary conditions were used, respectively. Free-stream turbulence intensity was set to 0.15%. All the simulations were performed at Cranfield's University Cluster (Delta). For each of the simulations, the HPC setup of 16 Ã 8 CPUS were used, and 12 CPU hours were needed.
An unstructured grid was used for the spatial discretisation of the domain. A refinement volume was defined in the vicinity of the wingtip vortex, to resolve the strong flow gradients within the vortex region (volume B in Figure 3) whose exact position was iteratively determined. Four grids were developed to examine the influence of the domain's discretization on the flow predictions with a total number of elements between 30.7•10 6 and 79.8•10 6 . The resolution across the region A was kept constant at 0.00833 C along all directions. Region B resolution ranged between 0.005833 C and 0.002333 C. (see Table 1). For the wing surface refinement, different spacings were examined on grid #3; namely a coarse (0.005833 C), a medium (0.004166 C) and a fine (0.0025 C) surface resolution.
Two turbulent models were used in this study; the Spalart-Allmaras (S.A.) and Reynold Stress Model (R.S.M.). As previously described by Pereira et al. 33 the R.S.M. approach shows accurate predictions for wingtip vortex flows, whereas the S.A., with the rotation correction of Dacles-Mariani et al., 43 was found to be the best compromise between accuracy and computational cost among a range of RANS closures based on the Boussinesq hypothesis. For the R.S.M., the sub-model of linearpressure stress model was used. For the near-wall boundary layer treatment, the standard wall function approach was adopted, respecting the condition of the first cell height, y þ > 30. For the simulations with the S.A. model, the turbulence production was modelled via the strain-vorticity-based production first introduced by Dacles-Mariani et al. 43 with a y + ≤1 in all near-wall regions.

Results and discussion CFD validation
The influence of region B spatial resolution, on the axial velocity component and static pressure coefficient at the vortex core, is shown in Figure 4. The wing's leading-edge was placed at x=c ¼ 0. The position of the trailing-edge at x=c ¼ 1 is indicated in Figure 4 by the dotted line. At x=c ¼ 1:66, the non-dimensional axial velocity at the vortex core was under-predicted with reference to the experimental data by 18.6%, 10.6% and 3.5% for grids #1, #2 and #3, respectively (Figure 4a). Grid #4 showed a discrepancy of around 4% from the test data up to x=c ¼ 1:23 after which the prediction was within 1% of the experiment. In terms of static pressure coefficient, a maximum discrepancy of 44%, 34%, 24% and 20% between the CFD and the test data were found at x=c ¼ 1:66, for grids #1, #2, #3 and #4, respectively (Figure 4b). Grids #3 and #4, showed a mean discrepancy of the vortex core static pressure coefficient of around 4% across the range of downstream positions. Craft et al. 32 and Pereira et al. 33 also reported the difficulty of modelling the static pressure coefficient at the vortex core in this test case. The influence of the grid's resolution along the spanwise direction (zaxis) on the vortex characteristics, is shown in Figure 5. Grids #3 and #4 were examined for this study as the finest and most accurate to predict the vortex core flow magnitudes in Figure 4. The comparison was performed at the first x=C ¼ 1:11, and last x=C ¼ 1:66, axial locations within the wake, for which experimental data were available. Between the two examined grids, the CFD model predicted the vortex structure independently of grid resolution. The agreement with the experimental results is also shown. Based on the above observations, grid #3 was chosen for further studies.
The static pressure coefficient C p , non-dimensional axial velocity u=U ∞ , and non-dimensional crossflow velocity U θ =U ∞ variations from grid #3, across the plane located at x=C ¼ 0:974, are shown in sub-figures Figure  6 Finally, the grid dependency of the static pressure distribution along the wing surface has been also studied, via the three different grid surface refinements of coarse, medium and fine, respectively. The results have been found to be independent of the surface grid selection and in agreement with the experimental data. The CFD returned near-field drag coefficients of 0.02954, 0.02935 and 0.02927 via the coarse, medium and fine grid resolution, respectively. These correspond to discrepancies of 0.9% and 0.3. For further analysis, the fine surface resolution was selected due to the negligible additional computational cost compared to coarse and medium grid resolutions.
A computational domain size study was also conducted to determine the influence of wind-tunnel walls onto the wingtip vortex characteristics and the appropriate domain size for the far-field, exergy-based drag analysis. The axes of y and z gradually increased until no effects were calculated (see also Figure 2). Namely, the wing's root remained mounted on the symmetry plane and every extension on y axis was applied both in upward and downward directions. The results of three domains are reported; domain #1 representing the baseline wind-tunnel size (see Figure  2), domain #2 in which z axis extended to 5 chords and   The analysis showed that the baseline domain size (domain #1, wind-tunnel representative) with x=C ¼ 7:5, y=C ¼ 0:66, and z=C ¼ 1 (see Figure 2), had a significant influence on the vortex strength. Increasing the computational domain size to domain #3 (x=C ¼ 7:5, y=C ¼ 20, and z=C ¼ 10) the non-dimensional axial velocity component at the vortex core reduced by 28% and the static pressure coefficient by 57% at x=C ¼ 1:11 (see Figure 7) in relation to domain #1. In addition, near-field drag reduced by 33% due to the reduced impact of the wingtip vortex on the wing. Between domains #2 (x=C ¼ 7:5, y=C ¼ 10, and z=C ¼ 5) and #3 no variations were found on the vortex strength, structure and wing's near-field drag. However, for the far-field, exergy-based analysis, the largest domain #3 was used to secure that all flow gradients, within the wake, have been integrated.
Spalart-Allmaras (S.A.) turbulence model was compared with the Reynolds Stress Model (R.S.M.) for domain #3. The non-dimensional axial velocity and static pressure coefficient at the vortex core are shown in Figure  8 where at x=C ¼ 1:66, R.S.M. compared to S.A., overpredicted the axial velocity, u=U ∞ and static pressure coefficient C p by 3.2% and 4.8%, respectively. Figure 8 also shows the variations between the two turbulent  models in the baseline computational domain #1. The impact of turbulence model on the non-dimensional crossflow velocity U θ =U ∞ is shown in Figure 9 where both turbulent models were found in good agreement in terms of vortex core location, strength and structure calculations. On the near-field drag calculations, a discrepancy of 1.5% between the two models was calculated. Namely, Spalart-Allmaras returned a drag coefficient value of 0.01934 and Reynolds Stress Model a value of 0.01964. Based on the above analysis, it is shown that Spalart-Allmaras model successfully predicted the wingtip vortex flow-field against the Reynolds Stress Model calculations in domain #3. In addition, following the restriction of the exergy balance method on the Boussinesq's hypothesis, for the far-field drag decomposition analysis, the turbulence model of Spalart-Allmaras, domain #3 and grid #3 have been adopted.
The comparison between the far-field, exergy-based and near-field drag predictions, as well as drag decomposition into exergy outflows and anergy terms (equation (3)), are shown in Figure 10. Each component of (equation (3)) was expressed in power counts (see (equation (12)), and then normalised by the near-field drag coefficient. The CFD calculated, near-field drag is shown by a solid line of constant value along the wake ( Figure  10). To enable fair comparisons, the near-field drag figure shown in Figure 10 was also produced using the domain #3 and grid #3 of the computational model. The discrepancy between the far-field and near-field total drag was 5.6% at the first wake point (x=C ¼ 1:1). One chord downstream, at x=C ¼ 2:1 this discrepancy was reduced to 3%, whereas at the last seven wake points (3:4≤x=C ≤ 4:1) was calculated between 1% and 0.55%.
Starting from the leading-edge of the wing, the exergy balance method provided information about the variation of drag and its components along the wing's solid body. At x=C ¼ 0:5, 68% of the total drag has been developed; 66% of which was linked with reversible phenomena and 34% with generated entropy (see Figure  10). It is indicated that a great amount of recovery potential has been induced from the interaction of the wing with the free-stream flow. Namely, in the wake it was calculated between 75% at x=C ¼ 1:1, and 68% at x=C ¼ 4:1 ( Figure 10). Entropy reflected approximately 30% of the total drag. The two components increased  almost linearly up to the 95% of the wing's chord (x=C ¼ 0:95), whereas a sharp decrease of 4% of the total exergy outflows was calculated at the trailing-edge ( Figure 10). This, however, was not reflected on an analogous increase in the anergy term (which increased only by 1%).
The theoretical formulation of exergy balance in equation (4), defines that the total far-field drag remains constant via the increase of dissipative terms of _ A Φ and _ A =T as a result of mechanical (_ ε m Þ and thermal (_ ε th Þ exergy outflows destruction. To quantify the numerical uncertainties that influence the total far-field drag calculation, equation (4) can be written as where dx denotes the variation of the axial position between two parallel transverse planes along the x axis. Since the total far-field drag must be independent on the spatial location of the transverse plane and the wave anergy has been deactivated, due to low-subsonic conditions, equation (13) is reduced to where _ E tot and _ A tot are the total exergy outflows and total anergy, respectively. However, in the case of finite number of cells, the right-hand side of equation is equal to the spatial variation of spurious drag, d dx Sp. Therefore, equation (14) yields: The terms of equation (15) can be expressed in power counts as explained in equation (12) and then divided by   the CFD calculated near-field drag coefficient. Thus, the spatial variation of spurious drag as a function of near-field drag can be shown ( Figure 11). Figure 11(a), shows the variations, within the wake, of the right-hand side components of equation (15), non-dimensionalised by the near-field drag. Higher rate of exergy destruction was calculated closer to the trailing-edge and reduced to the half on the next wake stations. The fluctuations on the results are due to the interpolation error linked with the 2D transverse plane locations. In contrary, the calculations on anergy production returned smoother variations due to the post-processing in volume integrals. Again, the higher rate was observed closer to the trailing-edge, whereas at the last computed point, (x=C ¼ 4:1), returned a magnitude of 0.05% of the total near-field drag.
The spatial variation of spurious drag nondimensionalised by the near-field drag is shown in Figure 11(b). The negative sign indicates the additional Figure 11. Evolution of exergy destruction and anergy production rates (a) and study of (in)consistency between the theoretical exergy balance formulation and CFD implementation (b). calculated exergy destruction against the total anergy terms increase. Due to the numerical fluctuations on the results, a trendline defined by a fourth order polynomial, demonstrates the variation of the spurious drag within the wake. In the first points, the difference between the exergy destruction and anergy increase, was around 0.4% of the total near-field drag. Two chords downstream, at x=C ¼ 3, it reached an average value of 0.15% and then reduced again to 0.05% at the last point (x=C ¼ 4:1). The higher spurious drag rate was calculated at the first wake stations due to the flow interaction with the blunt shape of the trailing-edge. A similar observation reported by Arntz et al. 26 where a transient region was shown, right downstream the aerodynamic body, where the flow 'reequilibrates itself' following the disturbance caused by the aerodynamic body. This also explained the asymptotic behaviour of the total far-field drag in Figure 10. Total drag decomposition into reversible and irreversible components via the exergy-based, far-field method, is shown in Figures 12 and 13 in terms of axial ð _ E u Þ and transverse ð _ E v Þ kinetic exergy, pressurework rate ð _ E p Þ and viscous dissipation ð _ A Φ Þ variations (as defined in equation (5) and equation (7)). The maximum magnitude of axial kinetic exergy ð _ E u Þ was found at the leading-edge ðx=C ¼ 0Þ and it was equivalent of 47% of the total near-field drag (left axis of Figure 12). At the same position, the pressure-work rate was found equivalent of 84% of the total near-field drag. Towards the wing's trailing-edge (x=C ¼ 1), due to pressure expansion, both _ E p and _ E u reduced almost linearly, up to the point of x=C ¼ 0:7, to 7.5% and 6%, respectively. The wingtip vortex development caused the increase of _ E u and _ E p magnitudes due to the static pressure drop and axial velocity increase within the vortex core. At the trailing-edge, the two terms calculated with the same magnitude and opposite sign. In the last examined wake station (x=C ¼ 4:1), due to dissipative phenomena, _ E u was found 0.7% and _ E p 0.05% of the near-field drag.
The variation of transverse kinetic exergy outflow _ E v , indicated the dominant impact of the lift-induced vortex on the wing's total drag (see in Figure 12 the right axis). The term represents the maximum amount of energy that  could be recuperated from the induced drag. At x=C ¼ 0:1, the transverse kinetic exergy represented the 30% of the total near-field drag. It increased almost linearly up to the trailing-edge, where the wingtip vortex was fully developed, returning a peak magnitude of 0:8D nf . The calculation highlighted the slow decay of the vortex within the wake. Namely, at the first wake station (x=C ¼ 1:1), _ E v reflected a magnitude equivalent of the 76% of the total drag and at the last wake point, x=C ¼ 4:1, the 68%. Due to low-subsonic conditions, thermal effects ( _ ϵ th ) had negligible impact in the total exergy-based drag (maximum magnitude at the trailingedge, C _ ϵ th ¼ 0.015 CD nf ), and therefore, the relative figure is not shown here.
The variation of viscous dissipation _ A Φ , within the wake of the wing is illustrated in Figure 13. The component increased almost linearly from the leading to the trailing-edge of the wing due to the boundary layer dissipation. At the leading-edge and due to the stagnation point, the viscous anergy was calculated 5% relative to the total near-field drag. In the wake, the dissipation increased with a slow rate due to the decay characteristics of the vortex. The calculated thermal anergy ( _ A =T ) was negligible and therefore is not shown here. In the far-field wake, beyond the demonstrated wake points, it is expected the total anergy to dominate the exergy balance terms (equation (4)) due to exergy outflows terms dissipation.
Flow-field analysis: Further considerations on wake recovery One of the major advantages of the far-field, exergy-based method, is the physical breakdown of the flow components responsible for the total drag of the wing and mapping of the sources associated with energy recovery. The mechanical exergy distributions are demonstrated in this section at a single crossflow plane within the wake as defined in Figure 14(a) aiming to quantify the potentially exploitable amount of wake energy by a closely integrated propulsion system. The transverse plane was located at the stream-wise distance, x=C ¼ 2:5 which is the location where the spatial rate of change of the drag components remains below 1% (see Figure 11a). The axial kinetic exergy outflow _ E u is shown in Figure 14(b), the transverse kinetic exergy _ E v in Figure 14(c) and the pressurework rate _ E p in Figure 14(d). In total, a potential propulsion system may recover, at this wake position (x=C ¼ 2:5), up to 2.8 kW of mechanical power predominantly linked with the transverse kinetic exergy terms (lift-induced drag) and at a much lesser extent with axial kinetic exergy (mainly induced by boundary layer) term. This outcome indicates that the wingtip is possibly a promising location lending itself to a closely coupled integration of a propulsion system designed to harvest part of the available exergy. The amount of maximum power recovery that could be achieved, depends on the design characteristics of the candidate propulsor. For instance, the maximum recovery for a pusher propulsor, with circular inlet, aligned with the axis of the vortex core and a diameter of 22% of the wing's span, would reflect to the 21% of the total amount of available power within the wake. Of course, other parameters such as installation effects, would significantly affect the flow features and therefore the system's recovery potential. However, the current analysis was focused on the applicability and understanding of method's potential benefits towards a more closely coupling between the aerodynamic body (wing and/or fuselage) and the propulsor.

Conclusions
In this work, the application of a far-field, exergy-based method for drag decomposition downstream of a subsonic 3D wing was demonstrated aiming to characterise the wake and determine the potentially recoverable amount of energy. Such an approach enables closer coupling between the wing and a propulsion system in order to recover part of the available wake energy which conventional propulsion integration methods don't usually provide. The study was based on an academic low aspect ratio wing, at an angle of attack of 10°under low-subsonic free-stream conditions of M = 0.15 whose induced flow-field is predominantly dominated by the formation of a strong wingtip vortex. This test enabled the development of the far-field method on a baseline, reference configuration before applying it to more representative geometries.
A RANS based CFD model was developed and validated against experimental data previously shown by Chow et al. 30 The validation was focused on the characteristics of the wingtip vortex generation and development within the wing's near-field but also across numerous axial positions along the wake. Prior to the application of the far-field drag decomposition method, a domain mesh and turbulence closure sensitivity studies were carried out. These studies determined the required size of the computational domain to mitigate the impact of the domain boundaries on the flow-field and also to ensure that the flow-field is fully developed within the computational domain. The required grid's spatial resolution to capture the flow gradients around the wingtip vortex was determined. It was finally shown that a Spallart-Allmaras turbulence closure offers a better compromise between accuracy and computational cost compared to a Reynolds Stress Model when domain size effects are mitigated.
Near-field drag calculations were used as reference to verify the accuracy of the far-field, exergy-based drag calculations. The far-field method was found to match the near-field total drag calculation within a maximum discrepancy of 5.6% at a distance of 0.1 wing's chords downstream the trailing-edge. Three chords downstream the trailing-edge exergy-based drag was over-calculated by 0.55%. Further drag decomposition along the wake was also carried out via the far-field method to determine the wake's power losses break down. The outcomes of the analysis showed that the lift-induced vortex together with the viscous dissipation, occurred mainly within the boundary layer, dominate the total drag. The methodology was able to provide a detailed physical drag decomposition and map the recoverable amount of energy included in the wake. It is thus shown that a far-field drag decomposition can potentially guide the design and integration of a propulsion system with the airframe from the very early stages of a development programme and support the specification of the key design characteristics of a propulsor to maximise the efficiency of the integrated system.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Engineering and Physical Sciences Research Council [I-CASE award number 16000148]; and Airbus.