# Energy-Based Aerodynamic Analysis on the Blended-Wing-Body Aircraft with Boundary Layer Ingestion

Authors: Wenyuan Zhao; Yanlai Zhang; Jianghao Wu
Journal: International Journal of Aerospace Engineering
 (2022)

---

## Abstract

This study is aimed at evaluating the energy saving benefit and its limitation of the boundary layer ingestion (BLI) applied to a transonic blended-wing-body (BWB) aircraft. The power balance method is adopted in the aerodynamic performance analysis of the aircraft. We further improve the one-dimensional analysis framework by adding the pressure volumetric work term in the power balance equation for the application in compressible flow with shock wave regions. A determining expression for the power-saving coefficient (PSC) of the BLI effect is also deduced, which is related to a BLI fraction, the airframe dissipation component ratios, and the specific thrust of the propulsor. This expression provides a way for the preliminary evaluation and comparison of BLI benefits applied in different aircraft configurations. The flow quantities of three different configurations are obtained by 3D Reynolds-averaged Navier–Stokes (RANS) numerical simulations for the calculation of aerodynamic dissipation components. The grid refinement and domain size study shows that the calculated force coefficients and total dissipation coefficient have a good convergence characteristic with the mesh refinement and are insensitive to the variation in domain size. The evaluated PSC of the BWB aircraft with a BLI propulsor is 5% using the previously developed determining expression. The maximum of PSC is 18.8%. The following comparison of dissipation components for the power-on and power-off cases shows that the benefit of BLI is slightly underestimated by the determining expression of PSC. The reason is that the variation in the boundary layer velocity profile at the trailing edge of the center body part is neglected, which is due to the installation of the propulsor. Furthermore, when the BLI propulsor is applied to a BWB aircraft, the variation in trailing edge vortex dissipation and volumetric pressure work term should be studied in detail.


---

## Body

## 1. Introduction

Several advanced subsonic transport aircraft configurations and related propulsion technologies have been investigated [1–6] to improve aircraft fuel efficiency and emission levels [7, 8]. Most of these new aircraft configurations integrate a propulsion system with an airframe by using an embedded installation method. In this manner, the propulsion system ingests a part of or all of the boundary layer/wake flow generated on the airframe. This airframe–propulsor integration technology is called boundary layer ingestion (BLI); it can potentially reduce total power consumption by up to 20% in preliminary studies [9]. Among the novel aircraft configurations that take advantage of BLI effects, blended-wing-body (BWB) aircraft, such as the Boeing BWB-450 [1], SAX-40 [2], and N3-X [3], feature wings merged with the fuselage and exhibit the advantage of higher airframe aerodynamic efficiency than conventional tube-and-wing aircraft [1]. The wide and flat upper surface of the BWB airframe enables the embedded installation of the propulsion system near the trailing edge, and a large part of the boundary layer flow generated on it can be ingested by the propulsion inlet. The estimated energy savings benefit of BWB aircraft with a BLI propulsor system is predicted to reach 3%–10% [10–13], which is affected by the fractions of boundary layer ingestion, the propulsor design parameters, and the application of active flow control equipment.Earlier studies of BLI concept applied to airplanes [14] used the estimation and reduction of airframe drag as the metric for quantifying the benefits of BLI. Subsequently, the estimations of the ram drag and pressure recovery, along with the drag of the airframe, were used as metrics for analyzing and optimizing the BLI configurations [10–13]. The aerodynamic and propulsive performance assessments in these earlier studies are based on a conventional thrust–drag bookkeeping method. In this method [15], the inner and outer flows are perfectly separated. A clear definition of thrust and drag of the power-on aircraft can be obtained by division of inner and outer flows. However, the inner flow cannot be separated from the outer flow in the BLI aircraft because the flow passing on the upper surface is also used as the intake flow of the propulsor. The clear definition of thrust and drag is compromised without the division of inner and outer flows. Thus, the aeropropulsive performance metric based on a thrust–drag bookkeeping is imprecise for the highly coupled BLI configurations.Due to the imprecision of the conventional thrust-drag bookkeeping method for the coupled configurations, Drela introduced the power balance method to analyze the flow field around these configurations [16]. In this method, the aerodynamic performance of the aircraft is analyzed by the balance between the mechanical energy provided by the propulsor and its consumption by the dissipation processes in the flow field around an aircraft. Next, total airframe dissipation is decomposed into physically intuitive components, which contribute to the understanding of the aeropropulsive interaction of the fluid. The total power requirement or total aerodynamic dissipation is selected as the aerodynamic performance metric in this method. Therefore, defining thrust or drag is no longer necessary, which is confusing in the BLI configurations.The power balance method has the potential to improve the understanding of the highly coupled aerodynamics of the BLI configurations. Analysis of the D8 aircraft based on this method has shown that the BLI configuration requires 6% to 9% less propulsive power relative to the non-BLI configuration under simulated cruise condition [17, 18]. The experiments done by Lv et al. [19, 20] found that the BLI propulsor uses the wake energy of body wake flow and reduces the wake energy of downstream wake flow, which decreases power consumption and simultaneously generates an equivalent net streamwise force compared to the non-BLI propulsor, under the same flight condition. The results also showed that wake ingestion and BLI configurations consume 10% and 18% less shaft power relative to the non-BLI propulsor, respectively [20]. Hall et al. summarized the BLI investigations on D8 and conducted a power balance-based control volume analysis of the BLI configuration [21]. They found that the primary cause of power saving from the BLI effects is reduction in wake dissipation and jet dissipation. The result of the control volume quantifying analysis is consistent with the results of the previously conducted computational fluid dynamics (CFD) simulation and experiment [21]. Uranga et al. emphasized the difference in dissipation of the boundary layer between the isolated configuration and the configuration with the BLI nacelle [22]. Sanders and Laskaridis identified the significant contributions of volumetric pressure work to the power balance of BLI configurations [23].The benefit of BLI varies with different configurations and propulsion systems. It is more evident when applied to a tube-and-wing airframe with small propulsors [6, 17, 18, 24] than when applied to a BWB airframe with distributed propulsion [12, 13]. In the BLI analysis of the D8 aircraft, the mechanical power requirement decreases when the mass flow rate of the propulsor increases under the same BLI fraction [21]. Furthermore, experimental studies on the bodies of revolution have found higher gains when boundary layer flow is completely ingested by the propulsor [25, 26]. These results show that the energy saving benefit of BLI varies due to different airframe dissipation characteristics, propulsion designs, and the BLI fraction. However, the relationship expression of this benefit with the three aspects has not yet been established. Through this expression, the design work of highly coupled configurations will be simplified, and a method for analyzing and comparing different configurations can be developed.The current study is aimed at evaluating the energy saving benefit and its limitation of BLI applied on a BWB aircraft by the power balance method. We further improve the one-dimensional analysis framework by adding the pressure volumetric work term in the power balance equation for the application in compressible flow with shock wave regions. We have deduced a determining expression for the energy saving benefit of BLI. Moreover, the expression will be related to the BLI fraction, airframe dissipation characteristics, and a propulsion design parameter. The aerodynamic dissipation components of the unpowered airframe were used in this expression. The limitation of the energy saving benefit of BLI can also be obtained by a specific set of the BLI fraction. This expression provides a way for the comparison of different aircraft configurations from the airframe and propulsor design perspectives. It is also appropriate for the analysis of different types of airframe-propulsor integration because it can capture the main effects of interaction.The remainder of the paper is structured as follows. In Section2, we simply present the power balance analysis method. Then, the determining expression of the power-saving coefficient (PSC) for BLI with the BLI fraction, the airframe dissipation characteristics, and the propulsor design parameter are presented to compare the BLI benefits applied to different aircraft. In Section 3, the flow simulation setup for unpowered and power-on aircraft are described. The detailed geometry description is first presented. Then, a propulsor model for obtaining the flow parameters at the inlet and outlet of the propulsor is constructed. Finally, the description of the integral volume for the calculation of dissipation components and the flow feature extraction method is presented. In Section 4, a study of the grid refinement, domain size, and integral volume size of the calculated dissipation and force coefficients is first performed. Then, the quantitative evaluation of PSC and its limitation of BLI applied on the BWB aircraft is presented, using the previously developed determining expression. The results of the BWB aircraft are also compared with those of another advanced aircraft (D8). Finally, a detailed dissipation component comparison between the BLI configuration and the non-BLI configuration of the BWB aircraft is performed to analyze the aeropropulsive coupling effects. The study findings are summarized in Section 5.


## 2. Energy-Based Aerodynamic Performance Analysis Method

In Section2.1, We describe the power balance relation in aerodynamic flow proposed by Drela [16]. That is, the mechanical power loss generated on the aircraft surface and in the wake/jet/vortex/wave regions must be balanced by the mechanical flow power provided by the propulsor. This theory method is suitable for analyzing aeropropulsive coupling configurations such as BLI aircraft because the airframe surface dissipation items do not evidently change with pressure field interference [16]. In Section 2.2, The total aerodynamic dissipation is decomposed into several components that are related to different flow features, which is helpful for the physical understanding of complex highly coupled aeropropulsive flows. In Section 2.3, some energy-based aerodynamic performance metrics are also introduced. In Section 2.4, We further improve the energy-based analysis framework firstly constructed by Hall et al. [21], through adding the pressure volumetric work term in the power balance equation for the application in compressible flow with shock wave regions. Finally, the determining expression of PSC with the BLI fraction, the airframe dissipation characteristics, and the propulsor design parameter are deduced for the primarily evaluation and comparison of different BLI aircraft.


### 2.1. Power Balance Relation in Aerodynamic Flow

We assume a steady flow in the current analysis. A cylinder control volume (CV) that includes a power-on aircraft is constructed for the analysis of the mechanical energy balance. The 2D slice of the control volume is shown in Figure1. The upstream and side parts of the outer boundary SO are placed sufficiently far from the aircraft to be under the free-stream condition. The side part is parallel to the free-stream velocity vector V⟶∞. The downstream part of SO, i.e., the Trefftz Plane (TP), is perpendicular to the free-stream velocity vector V⟶∞. The inner boundary SB covers the body surface with cut planes at the propulsor inlet and outlet locations to exclude the inflow of the ducted propulsor from the control volume. The inner boundary SB has no moving part. Any small mass added by the propulsion system is ignored. As derived by Drela [16], the mechanical power balance equation can be described as
(1)PV+PK=Φ+Ė−FXV∞.Figure 1
2D slice of the 3D control volume used for the power balance analysis.The two terms on the left side of Equation (1) are the production and inflow. The three terms on the right side are the mechanical power consumption caused by viscosity, irreversible outflow, and the power consumption used to increase the aircraft’s potential energy.The first termPV on the left side of Equation (1) is the net pressure-volume power, provided by the fluid expanding against ambient pressure.
(2)PV≡∭CVp−p∞∇⋅VdV.This term is equivalent to zero in incompressible flow, while the contribution of this term is significant in compressible flow [23]. It is usually assumed to be reversible, except in regions containing shocks.The second termPK on the left side of Equation (1) is the net mechanical power into CV. It is equal to the flow rate of mechanical energy through the inner boundary SB, which, in turn, is equal to zero except in the inlet and outlet planes of the propulsor.
(3)PK≡∯SB−p−p∞+12ρV2−V∞2V⋅n^dSB.It is the sum of net pressure work rate and kinetic energy flow rate across the inner boundary into the CV.PK can be related to fan shaft power and fuel consumption. If information about engine efficiency cannot be obtained, then it may function as a metric for the energy or fuel consumption.The third term−FXV∞ on the right side of Equation (1) is the mechanical power consumed or provided during the climb or glide of the aircraft, which is equal to Wḣ (the weight of the aircraft W multiplied by the climb rate ḣ) due to the relation of the weight of the aircraft during the steady climb or glide without banks. FX denotes the net streamwise aerodynamic force on the power-on aircraft and is equivalent to the inner boundary net streamwise force and momentum flow FBX or outer-boundary streamwise force and momentum flow FOX due to the integral momentum equation.
(4)FX=FBX=FOX,FBX=∯SBpn^−τ+ρVV⋅n^⋅exdSB,FOX=∯SO−p−p∞n^+ρV−V∞V⋅n^⋅exdSO.Thex-axis is along the free-stream velocity vector V∞. When the aircraft cruises at stable altitude and velocity, FX=0. −FXV∞ is a reversible energy outflow.The first termΦ on the right side of Equation (1) quantifies the rate at which mechanical energy is consumed by viscous dissipation within the control volume.
(5)Φ≡∭CVτ¯¯⋅∇⋅VdV.The second termĖ on the right side of Equation (1) is the irreversible loss of mechanical energy from the outflow due to flow nonuniformity on the TP.
(6)Ė≡∬TP12ρu2+v2+w2V∞+u+p−p∞udSOTP.When the TP is moved further downstream, this part of the kinetic energy is converted into thermal energy through viscous dissipation within the shear layer and turbulent mixing. Then, the total aerodynamic dissipationΦ∗ can be defined as
(7)Φ∗≡Φ+Ė−PV,which is built up by dissipation in different flow features such as shear layer, trailing edge vortex propulsive jet, and shock wave.Equation (1) can be rearranged as
(8)−FXV∞=PK−Φ∗,which is similar to the traditional momentum based thrust–drag bookkeeping system. In the latter system, net streamwise force of a power-on aircraftFX is divided into drag and thrust parts. Then, the two components are estimated in separated flows called external flow and inner flow. The interaction effect of external and inner flow on thrust–drag bookkeeping is neglected in the momentum-based system, which is inappropriate for highly coupled aeropropulsive flow. In the energy-based analysis system, the decomposition of net streamwise force is translated to the decomposition of net mechanical power residual −FXV∞. Given PK and FX of a power-on aircraft, the total aerodynamic dissipation Φ∗ can be extracted from Equation (8). The interaction effect of external flow and inner flow is considered by calculating aerodynamic dissipation components and power supplement component in the same coupled aeropropulsive flow, which makes the energy-based method more appropriate for the aerodynamic analysis of the highly coupled aeropropulsive flow.


### 2.2. Aerodynamic Dissipation Breakdown

The second term on the right-hand side of Equation (7) can be divided into three components:
(9)Ė=Ėa+Ėv+Ėp,each of which will be introduced in detail next. The result is a physical insightful decomposition of the total aerodynamic dissipationΦ∗:
(10)Φ∗≡Φ+Ėa+Ėv+Ėp−PV,where the first term on the right-hand side of Equation (10) represents viscous dissipation that occurred in the control volume described in Section 2.1, the next three terms are the three irreversible energies flowing out through TP components, and the last term is pressure work related to isotropic deforming of the fluid, which is only obvious when there is a shock wave in the flow field. Heating and thermal conduction are not considered in this paper.Wake streamwise kinetic energy deposition rateĖa is as follows:
(11)Ėa=∬TP12ρu2V∞+udSOTP.This represents the streamwise kinetic energy depositing rate of the TP outflow. It is always positive. This term decays quickly via the dissipation progress in the wake and jet.Wake streamwise kinetic energy deposition rateĖv is as follows:
(12)Ėv=∬TP12ρv2+w2V∞+udSOTP.It is the transverse kinetic energy depositing rate of the TP outflow, which decays with the dissipation generated in vortex mixing.Wake pressure-defect work rateĖp is as follows:
(13)Ėp=∬TPp−p∞udSOTP.It is the rate of pressure work done on the fluid crossing the TP at some pressurep different from the freestream pressure p∞. This item decays to a small amount at some TP with a distance of about one or two character-chord-lengths from the trailing edge, where the static pressure recovers to ambient pressure p∞.Moving the TP from the trailing edge of the aircraft to the far downstream, theĖa, Ėv, and Ėp components decay and the mechanical power is transferred from them to the inner energy of the fluid through the viscous dissipation Φ and fluid volumetric pressure work PV in the added control volume, as shown in Figure 2. By placing the TP within a small distance from the trailing edge of the aircraft, potential-near-field transients are generated in the Ėa, Ėv, and Ėp components, which can be canceled in the overall Ė sum [16]. By choosing a TP position where the static pressure is recovered to ambient pressure p∞, the potential-near-field-related parts in the individual Ėa, Ėv, and Ėp components disappear, and the remaining parts can be related to different flow features in the downstream flow.Figure 2
Distribution of different dissipation sources.By placing the TP sufficiently far from the aircraft, all dissipation progress have been completed. For a powered-on aircraft, the total aerodynamic dissipation can be decomposed as(14)Φ∗=Φbl+Φwake+Φvortex+Φjet−PV,where the dissipation components on the right side account for the dissipations generated in the boundary layers of the body surface, wake, trailing edge vortex, and propulsion system jet, as labeled in Figure2. −PV is the pressure work related to the isotropic deformation of the fluid, which is only irreversible in the shock wave region [23].For an unpowered aircraft, Equation (14) is reduced to
(15)Φaf∗′=Φbl′+Φwake′+Φvortex′−PV′.The profile dissipationΦp′ in an unpowered aircraft is defined as
(16)Φp′=Φ′bl+Φwake′,where prime ()′ refers to unpowered aircraft characteristics. Φp′ contains dissipation in the surface boundary layers and wakes. Notably, Φp′ is calculated with the TP located at a downstream position where static pressure is recovered to ambient static pressure p∞.When the static pressure is recovered to the ambient pressurep∞, and the streamwise perturbation velocity approximates zero, the mechanical energy deposited in the TP at this streamwise location (p=p∞) is actually the contribution of the transverse perturbation velocity. Then, it is gradually converted to inner energy by the vortex dissipation [16]. Thus, the vortex dissipation Φ′vortex is equivalent to the mechanical energy deposited in the TP at this streamwise location (p=p∞),
(17)Φ′vortex=Ėv.


### 2.3. Aerodynamic Performance Metrics Defined by Energy-Based Parameters

The lift-to-drag ratioK is usually chosen as the aerodynamic efficiency metric for the aircraft. However, the force-based aircraft drag is difficult to obtain when complex aeropropulsive coupling effects are exerted around aircraft with a highly integrated propulsor. Thus, the lift-to-drag ratio K is also inaccurate. A similar aerodynamic efficient parameter should be proposed based on power balance analysis. We call it the lift-power-to-airframe-dissipation ratio, and the definition expression is
(18)KE=LV∞Φ∗−Φjet,withΦjet representing jet dissipation.In this part, PSC is chosen as the metric for the performance comparison of BLI configuration with non-BLI configuration. With regard to the reduction in power consumption caused by BLI, the PSC was first introduced by Smith [9]. This metric is suitable for comparing the aerodynamic performance of different airframe–propulsor-integrated configurations because it has a consistent energy point of view with the power balance method. Since most of the fuel is consumed during the cruise throughout the flight envelope, the PSC is usually calculated in the cruise condition (FX=0). PSC is defined as
(19)PSC≡PK′−PKPK′.


### 2.4. Power Balance Analysis Framework of the BLI Aircraft

A simplified power balance analysis framework should be deduced to provide the power requirement constraint of the propulsor model described in Section3.2. for obtaining the flow conditions at the propulsor inlet and outlet used in the power-on flow simulation of the BWB aircraft with a BLI propulsor. As is proven by Hall et al. [21], the shear layer dissipation components are insensitive to static pressure changes brought by the ingesting propulsor. Thus, the unpowered airframe dissipation components can be used to assess the performance of the BLI configuration.With the assumption of a uniform jet speedVjet at the free-stream static pressure just downstream of the propulsor outlet, the mechanical flow power PK provided by the propulsor with BLI effects can be expressed as [21]
(20)PK=12Vjet2−V∞2mf+Φ′inl,where the first term on the right-hand side is the mechanical energy flux at the propulsor outlet and the second term is that at the propulsor inlet, equal to the dissipationΦ′inl that occurs upstream of the inlet within the ingested streamtube.The jet dissipation can be simplified as(21)Φjet=12Vjet−V∞2mf.The total aerodynamic dissipationΦ∗ in Equation (8) can be divided into airframe dissipation and jet dissipation:
(22)Φ∗=Φaf∗+Φjet.With the knowledge that the energy-saving benefit of BLI is the reduction inΦwake and Φjet [21] and the assumptions that the dissipation generated on the aircraft surface does not change evidently with the propulsion system and the vortex dissipation does not change due to the installation of the propulsor at a given lift coefficient [21], by introducing the BLI fraction fBLI and the wake dissipation ratio α, the total dissipation of the BLI aircraft can be estimated in terms of the dissipation components of the unpowered airframe:
(23)Φaf∗=Φaf∗′−fBLIΦwake′=1−fBLIαΦaf∗′,wherefBLIΦwake′ represents the reduction of wake dissipation due to BLI and α represents the ratio of wake dissipation to total airframe dissipation of the unpowered aircraft. Here, Hall et al. [21] assumed that the ratio of the reduction of wake dissipation to the total wake dissipation of the unpowered aircraft is equal to the BLI fraction fBLI.The BLI fractionfBLI is defined as
(24)fBLI=Φ′inlΦbl′,which means the ratio of the dissipation that occurs upstream of the inlet within the ingested streamtubeΦ′inl to that that occurs in all of the flow from the airframe surface boundary layer Φ′bl.The wake dissipation ratioα is defined as
(25)α=Φwake′Φaf∗′.By combining Equations (20), (21), (22), (23), (24), and (25), the power balance equation, Equation (8) at FX=0, can be rewritten as
(26)mfV∞2VjetV∞−1−1−fBLIβΦaf∗′=0,which is the energy requirement constraint of the propulsor model. The parameter β represents the ratio of the profile dissipation to total airframe dissipation, which is defined as
(27)β=Φp′Φaf∗′.Φp′ is the profile dissipation of the unpowered airframe, which consists of dissipation in the surface boundary layer and trailing edge wake.By introducing Equation (26) into Equation (21), the jet dissipation can be rewritten as
(28)Φjet=12ε1−fBLIβ2Φaf∗′,where the propulsor design parameter ε is defined as
(29)ε=Φaf∗′mfV∞2.For the unpowered airframe,Φaf∗′ equals to the conventional drag power DV∞. In the cruise condition (FX=0), the drag power DV∞ is equal to the thrust power TV∞, which can be expressed as mfV∞Vjet−V∞. Thus, ε is equivalent to the conventional specific thrust Vjet/V∞−1.By lettingfBLI=0, the jet dissipation of non-BLI aircraft Φjet′ can be rewritten as
(30)Φjet′=12εΦaf∗′.By combing Equations (8), (22), (23), and (28), the total energy requirement PK in FX=0 can be expressed as
(31)PK=1−fBLIα+12ε1−fBLIβ2Φaf∗′.By lettingfBLI=0, the total energy requirement of the non-BLI configuration can be expressed as
(32)PK′=1+12εΦaf∗′.With the assumption of an equivalent mass flow rateṁ between the BLI configuration and the non-BLI configuration, the PSC defined by Equation (19) brought about by BLI can be expressed as
(33)PSC=fBLIα+1/2ε2fBLIβ−fBLIβ21+1/2ε.PSC is determined on the basis of four factors: the BLI fractionfBLI, the wake dissipation ratioα, the profile dissipation ratio β, and the propulsion design parameter ε. With a higher BLI fraction fBLI and a larger ε, PSC will be evident and give known values of α and β. Given the value of α=9% and β=60% of the BWB aircraft designed by Yan et al. [27], the variation of PSC with fBLI and ε is shown in Figure 3.Figure 3
The variation of PSC withfBLI and ε.SettingfBLI=1, the maximum of PSC can be expressed as
(34)PSCmaximum=α+1/2ε2β−β21+1/2ε.The limitation of BLI benefits at a givenṁ is determined on the basis of the wake dissipation ratio α, the profile dissipation ratio β, and the propulsion specified thrust ε. Given similar values of α and β, PSCmaximum is largely affected by the design parameter of the propulsion system ε. A larger PSCmaximum can be achieved with a higher ε.According to Equation (33), the PSC consists of two terms: one is for the reduction of aircraft dissipation and the other is from jet dissipation reduction.Combining Equations (28), (30), and (32), the term from jet dissipation reduction PSCjet is
(35)PSCjet=ΔΦjetPK′=Φjet′−ΦjetPK′=1/2ε1+1/2ε2fBLIβ−fBLIβ2,which is a function ofε, fBLI, andβ. The value of PSCjet increases monotonically with ε and fBLIβ, since fBLI≤1 and β≤1.Combining Equations (23) and (32), the term from airframe dissipation reduction PSCaf is
(36)PSCaf=ΔΦafPK′=Φaf′−ΦafPK′=11+1/2εfBLIα,which is the function ofε, fBLI, and α. The value of PSCaf increases with fBLIα but decreases with ε.


## 2.1. Power Balance Relation in Aerodynamic Flow

We assume a steady flow in the current analysis. A cylinder control volume (CV) that includes a power-on aircraft is constructed for the analysis of the mechanical energy balance. The 2D slice of the control volume is shown in Figure1. The upstream and side parts of the outer boundary SO are placed sufficiently far from the aircraft to be under the free-stream condition. The side part is parallel to the free-stream velocity vector V⟶∞. The downstream part of SO, i.e., the Trefftz Plane (TP), is perpendicular to the free-stream velocity vector V⟶∞. The inner boundary SB covers the body surface with cut planes at the propulsor inlet and outlet locations to exclude the inflow of the ducted propulsor from the control volume. The inner boundary SB has no moving part. Any small mass added by the propulsion system is ignored. As derived by Drela [16], the mechanical power balance equation can be described as
(1)PV+PK=Φ+Ė−FXV∞.Figure 1
2D slice of the 3D control volume used for the power balance analysis.The two terms on the left side of Equation (1) are the production and inflow. The three terms on the right side are the mechanical power consumption caused by viscosity, irreversible outflow, and the power consumption used to increase the aircraft’s potential energy.The first termPV on the left side of Equation (1) is the net pressure-volume power, provided by the fluid expanding against ambient pressure.
(2)PV≡∭CVp−p∞∇⋅VdV.This term is equivalent to zero in incompressible flow, while the contribution of this term is significant in compressible flow [23]. It is usually assumed to be reversible, except in regions containing shocks.The second termPK on the left side of Equation (1) is the net mechanical power into CV. It is equal to the flow rate of mechanical energy through the inner boundary SB, which, in turn, is equal to zero except in the inlet and outlet planes of the propulsor.
(3)PK≡∯SB−p−p∞+12ρV2−V∞2V⋅n^dSB.It is the sum of net pressure work rate and kinetic energy flow rate across the inner boundary into the CV.PK can be related to fan shaft power and fuel consumption. If information about engine efficiency cannot be obtained, then it may function as a metric for the energy or fuel consumption.The third term−FXV∞ on the right side of Equation (1) is the mechanical power consumed or provided during the climb or glide of the aircraft, which is equal to Wḣ (the weight of the aircraft W multiplied by the climb rate ḣ) due to the relation of the weight of the aircraft during the steady climb or glide without banks. FX denotes the net streamwise aerodynamic force on the power-on aircraft and is equivalent to the inner boundary net streamwise force and momentum flow FBX or outer-boundary streamwise force and momentum flow FOX due to the integral momentum equation.
(4)FX=FBX=FOX,FBX=∯SBpn^−τ+ρVV⋅n^⋅exdSB,FOX=∯SO−p−p∞n^+ρV−V∞V⋅n^⋅exdSO.Thex-axis is along the free-stream velocity vector V∞. When the aircraft cruises at stable altitude and velocity, FX=0. −FXV∞ is a reversible energy outflow.The first termΦ on the right side of Equation (1) quantifies the rate at which mechanical energy is consumed by viscous dissipation within the control volume.
(5)Φ≡∭CVτ¯¯⋅∇⋅VdV.The second termĖ on the right side of Equation (1) is the irreversible loss of mechanical energy from the outflow due to flow nonuniformity on the TP.
(6)Ė≡∬TP12ρu2+v2+w2V∞+u+p−p∞udSOTP.When the TP is moved further downstream, this part of the kinetic energy is converted into thermal energy through viscous dissipation within the shear layer and turbulent mixing. Then, the total aerodynamic dissipationΦ∗ can be defined as
(7)Φ∗≡Φ+Ė−PV,which is built up by dissipation in different flow features such as shear layer, trailing edge vortex propulsive jet, and shock wave.Equation (1) can be rearranged as
(8)−FXV∞=PK−Φ∗,which is similar to the traditional momentum based thrust–drag bookkeeping system. In the latter system, net streamwise force of a power-on aircraftFX is divided into drag and thrust parts. Then, the two components are estimated in separated flows called external flow and inner flow. The interaction effect of external and inner flow on thrust–drag bookkeeping is neglected in the momentum-based system, which is inappropriate for highly coupled aeropropulsive flow. In the energy-based analysis system, the decomposition of net streamwise force is translated to the decomposition of net mechanical power residual −FXV∞. Given PK and FX of a power-on aircraft, the total aerodynamic dissipation Φ∗ can be extracted from Equation (8). The interaction effect of external flow and inner flow is considered by calculating aerodynamic dissipation components and power supplement component in the same coupled aeropropulsive flow, which makes the energy-based method more appropriate for the aerodynamic analysis of the highly coupled aeropropulsive flow.


## 2.2. Aerodynamic Dissipation Breakdown

The second term on the right-hand side of Equation (7) can be divided into three components:
(9)Ė=Ėa+Ėv+Ėp,each of which will be introduced in detail next. The result is a physical insightful decomposition of the total aerodynamic dissipationΦ∗:
(10)Φ∗≡Φ+Ėa+Ėv+Ėp−PV,where the first term on the right-hand side of Equation (10) represents viscous dissipation that occurred in the control volume described in Section 2.1, the next three terms are the three irreversible energies flowing out through TP components, and the last term is pressure work related to isotropic deforming of the fluid, which is only obvious when there is a shock wave in the flow field. Heating and thermal conduction are not considered in this paper.Wake streamwise kinetic energy deposition rateĖa is as follows:
(11)Ėa=∬TP12ρu2V∞+udSOTP.This represents the streamwise kinetic energy depositing rate of the TP outflow. It is always positive. This term decays quickly via the dissipation progress in the wake and jet.Wake streamwise kinetic energy deposition rateĖv is as follows:
(12)Ėv=∬TP12ρv2+w2V∞+udSOTP.It is the transverse kinetic energy depositing rate of the TP outflow, which decays with the dissipation generated in vortex mixing.Wake pressure-defect work rateĖp is as follows:
(13)Ėp=∬TPp−p∞udSOTP.It is the rate of pressure work done on the fluid crossing the TP at some pressurep different from the freestream pressure p∞. This item decays to a small amount at some TP with a distance of about one or two character-chord-lengths from the trailing edge, where the static pressure recovers to ambient pressure p∞.Moving the TP from the trailing edge of the aircraft to the far downstream, theĖa, Ėv, and Ėp components decay and the mechanical power is transferred from them to the inner energy of the fluid through the viscous dissipation Φ and fluid volumetric pressure work PV in the added control volume, as shown in Figure 2. By placing the TP within a small distance from the trailing edge of the aircraft, potential-near-field transients are generated in the Ėa, Ėv, and Ėp components, which can be canceled in the overall Ė sum [16]. By choosing a TP position where the static pressure is recovered to ambient pressure p∞, the potential-near-field-related parts in the individual Ėa, Ėv, and Ėp components disappear, and the remaining parts can be related to different flow features in the downstream flow.Figure 2
Distribution of different dissipation sources.By placing the TP sufficiently far from the aircraft, all dissipation progress have been completed. For a powered-on aircraft, the total aerodynamic dissipation can be decomposed as(14)Φ∗=Φbl+Φwake+Φvortex+Φjet−PV,where the dissipation components on the right side account for the dissipations generated in the boundary layers of the body surface, wake, trailing edge vortex, and propulsion system jet, as labeled in Figure2. −PV is the pressure work related to the isotropic deformation of the fluid, which is only irreversible in the shock wave region [23].For an unpowered aircraft, Equation (14) is reduced to
(15)Φaf∗′=Φbl′+Φwake′+Φvortex′−PV′.The profile dissipationΦp′ in an unpowered aircraft is defined as
(16)Φp′=Φ′bl+Φwake′,where prime ()′ refers to unpowered aircraft characteristics. Φp′ contains dissipation in the surface boundary layers and wakes. Notably, Φp′ is calculated with the TP located at a downstream position where static pressure is recovered to ambient static pressure p∞.When the static pressure is recovered to the ambient pressurep∞, and the streamwise perturbation velocity approximates zero, the mechanical energy deposited in the TP at this streamwise location (p=p∞) is actually the contribution of the transverse perturbation velocity. Then, it is gradually converted to inner energy by the vortex dissipation [16]. Thus, the vortex dissipation Φ′vortex is equivalent to the mechanical energy deposited in the TP at this streamwise location (p=p∞),
(17)Φ′vortex=Ėv.


## 2.3. Aerodynamic Performance Metrics Defined by Energy-Based Parameters

The lift-to-drag ratioK is usually chosen as the aerodynamic efficiency metric for the aircraft. However, the force-based aircraft drag is difficult to obtain when complex aeropropulsive coupling effects are exerted around aircraft with a highly integrated propulsor. Thus, the lift-to-drag ratio K is also inaccurate. A similar aerodynamic efficient parameter should be proposed based on power balance analysis. We call it the lift-power-to-airframe-dissipation ratio, and the definition expression is
(18)KE=LV∞Φ∗−Φjet,withΦjet representing jet dissipation.In this part, PSC is chosen as the metric for the performance comparison of BLI configuration with non-BLI configuration. With regard to the reduction in power consumption caused by BLI, the PSC was first introduced by Smith [9]. This metric is suitable for comparing the aerodynamic performance of different airframe–propulsor-integrated configurations because it has a consistent energy point of view with the power balance method. Since most of the fuel is consumed during the cruise throughout the flight envelope, the PSC is usually calculated in the cruise condition (FX=0). PSC is defined as
(19)PSC≡PK′−PKPK′.


## 2.4. Power Balance Analysis Framework of the BLI Aircraft

A simplified power balance analysis framework should be deduced to provide the power requirement constraint of the propulsor model described in Section3.2. for obtaining the flow conditions at the propulsor inlet and outlet used in the power-on flow simulation of the BWB aircraft with a BLI propulsor. As is proven by Hall et al. [21], the shear layer dissipation components are insensitive to static pressure changes brought by the ingesting propulsor. Thus, the unpowered airframe dissipation components can be used to assess the performance of the BLI configuration.With the assumption of a uniform jet speedVjet at the free-stream static pressure just downstream of the propulsor outlet, the mechanical flow power PK provided by the propulsor with BLI effects can be expressed as [21]
(20)PK=12Vjet2−V∞2mf+Φ′inl,where the first term on the right-hand side is the mechanical energy flux at the propulsor outlet and the second term is that at the propulsor inlet, equal to the dissipationΦ′inl that occurs upstream of the inlet within the ingested streamtube.The jet dissipation can be simplified as(21)Φjet=12Vjet−V∞2mf.The total aerodynamic dissipationΦ∗ in Equation (8) can be divided into airframe dissipation and jet dissipation:
(22)Φ∗=Φaf∗+Φjet.With the knowledge that the energy-saving benefit of BLI is the reduction inΦwake and Φjet [21] and the assumptions that the dissipation generated on the aircraft surface does not change evidently with the propulsion system and the vortex dissipation does not change due to the installation of the propulsor at a given lift coefficient [21], by introducing the BLI fraction fBLI and the wake dissipation ratio α, the total dissipation of the BLI aircraft can be estimated in terms of the dissipation components of the unpowered airframe:
(23)Φaf∗=Φaf∗′−fBLIΦwake′=1−fBLIαΦaf∗′,wherefBLIΦwake′ represents the reduction of wake dissipation due to BLI and α represents the ratio of wake dissipation to total airframe dissipation of the unpowered aircraft. Here, Hall et al. [21] assumed that the ratio of the reduction of wake dissipation to the total wake dissipation of the unpowered aircraft is equal to the BLI fraction fBLI.The BLI fractionfBLI is defined as
(24)fBLI=Φ′inlΦbl′,which means the ratio of the dissipation that occurs upstream of the inlet within the ingested streamtubeΦ′inl to that that occurs in all of the flow from the airframe surface boundary layer Φ′bl.The wake dissipation ratioα is defined as
(25)α=Φwake′Φaf∗′.By combining Equations (20), (21), (22), (23), (24), and (25), the power balance equation, Equation (8) at FX=0, can be rewritten as
(26)mfV∞2VjetV∞−1−1−fBLIβΦaf∗′=0,which is the energy requirement constraint of the propulsor model. The parameter β represents the ratio of the profile dissipation to total airframe dissipation, which is defined as
(27)β=Φp′Φaf∗′.Φp′ is the profile dissipation of the unpowered airframe, which consists of dissipation in the surface boundary layer and trailing edge wake.By introducing Equation (26) into Equation (21), the jet dissipation can be rewritten as
(28)Φjet=12ε1−fBLIβ2Φaf∗′,where the propulsor design parameter ε is defined as
(29)ε=Φaf∗′mfV∞2.For the unpowered airframe,Φaf∗′ equals to the conventional drag power DV∞. In the cruise condition (FX=0), the drag power DV∞ is equal to the thrust power TV∞, which can be expressed as mfV∞Vjet−V∞. Thus, ε is equivalent to the conventional specific thrust Vjet/V∞−1.By lettingfBLI=0, the jet dissipation of non-BLI aircraft Φjet′ can be rewritten as
(30)Φjet′=12εΦaf∗′.By combing Equations (8), (22), (23), and (28), the total energy requirement PK in FX=0 can be expressed as
(31)PK=1−fBLIα+12ε1−fBLIβ2Φaf∗′.By lettingfBLI=0, the total energy requirement of the non-BLI configuration can be expressed as
(32)PK′=1+12εΦaf∗′.With the assumption of an equivalent mass flow rateṁ between the BLI configuration and the non-BLI configuration, the PSC defined by Equation (19) brought about by BLI can be expressed as
(33)PSC=fBLIα+1/2ε2fBLIβ−fBLIβ21+1/2ε.PSC is determined on the basis of four factors: the BLI fractionfBLI, the wake dissipation ratioα, the profile dissipation ratio β, and the propulsion design parameter ε. With a higher BLI fraction fBLI and a larger ε, PSC will be evident and give known values of α and β. Given the value of α=9% and β=60% of the BWB aircraft designed by Yan et al. [27], the variation of PSC with fBLI and ε is shown in Figure 3.Figure 3
The variation of PSC withfBLI and ε.SettingfBLI=1, the maximum of PSC can be expressed as
(34)PSCmaximum=α+1/2ε2β−β21+1/2ε.The limitation of BLI benefits at a givenṁ is determined on the basis of the wake dissipation ratio α, the profile dissipation ratio β, and the propulsion specified thrust ε. Given similar values of α and β, PSCmaximum is largely affected by the design parameter of the propulsion system ε. A larger PSCmaximum can be achieved with a higher ε.According to Equation (33), the PSC consists of two terms: one is for the reduction of aircraft dissipation and the other is from jet dissipation reduction.Combining Equations (28), (30), and (32), the term from jet dissipation reduction PSCjet is
(35)PSCjet=ΔΦjetPK′=Φjet′−ΦjetPK′=1/2ε1+1/2ε2fBLIβ−fBLIβ2,which is a function ofε, fBLI, andβ. The value of PSCjet increases monotonically with ε and fBLIβ, since fBLI≤1 and β≤1.Combining Equations (23) and (32), the term from airframe dissipation reduction PSCaf is
(36)PSCaf=ΔΦafPK′=Φaf′−ΦafPK′=11+1/2εfBLIα,which is the function ofε, fBLI, and α. The value of PSCaf increases with fBLIα but decreases with ε.


## 3. Configuration Description and Numerical Simulation

To analyze the aerodynamic dissipation constitution of the BWB aircraft with BLI propulsors, 3D Reynolds-averaged Navier–Stokes (RANS) numerical simulations are performed for three different configurations, namely,C0, C1, and C2. Detailed descriptions of the three aircraft geometries are presented in the following subsection. Flow simulations are conducted using a second-order finite volume solver, and the Spalart–Allmaras turbulence model is selected. An average wall y+ of less than 1 is used. Half of the configuration is modeled with a symmetry plane along the centerline. As is shown in Figure 4, the pressure-far-field boundary condition is used at the outer boundary of the domain, which prescribes the static conditions of the free stream, the Mach number, and the direction of velocity. The propulsor inlet plane is modeled with the pressure-outlet boundary condition with the specified target mass flow rate and velocity direction. Meanwhile, the exit plane is modeled with the mass-flow-inlet boundary condition with the specified stagnation temperature, mass flow rate, and velocity direction.Figure 4
Boundary conditions used in the flow simulation.The detailed configuration description is presented in Section3.1. A propulsor model is constructed in Section 3.2 to obtain the flow parameters at the inlet and outlet of the propulsor. The description of the integral volume for the calculation of dissipation components and the flow feature extraction method is presented in Section 3.3.


### 3.1. Aircraft Configuration Description

We consider a subsonic BWB transport aircraft designed by Yan et al. [27], which has a 6,000 nm air range with a cruise Mach of 0.85 and a cruise altitude of 35,000 ft. It can carry 350 passengers with a 507,495 lb maximum takeoff gross weight. The designed lift coefficient in the cruise condition is CL=0.36. Three different configurations are constructed, namely, C0, C1, and C2, to analyze the effects of aeropropulsor coupling on the aerodynamic performance of the BWB aircraft. C0 is a clean BWB airframe with no propulsor installed on it, which is actually the unpowered configuration. The BWB airframe is composed of three parts: center body, blend wing, and outer wing. C1 consists of a BWB airframe equivalent to C0 and a propulsor placed over the rear part of the upper surface of the center body of the airframe near the trailing edge. C1 represents the power-on configuration with a non-BLI propulsor. C2 has the same BWB airframe as C0 and C1 and a BLI propulsor embedded in the rear part of the upper surface of the center body of the airframe to ingest the boundary layer flow generated on the airframe surface. The detailed vertical and symmetry plane views of the three configurations are shown in Figure 5.Figure 5
Configurations used in the energy-based aerodynamic performance analysis.3D RANS calculations are performed for the three configurations. The free-stream Mach number is set to 0.85; and the Reynolds number is approximately7×107, which is the same under the aircraft cruise condition. Given that the BLI propulsion system is embedded in the rear part of the upper surface of the center body, the lift coefficients CL of C1 and C2 may differ from that of C0 at the same attack angle. Thus, we perform the calculation on a range of attack angles and compare the dissipation components of the three configurations with a similar lift coefficient and not the same angle. The initial values of flow parameters at the propulsor inlet and outlet of C2 are calculated using the component-based propulsor model described in Section 3.2. The aircraft dissipation components used in the component-based propulsor model are the results of the unpowered configuration C0. Then, a calculation matrix constructed by the attack angle array and the fan operation array is generated, which is used to find the cruise point (CL=0.36,CX=0) of C2. The flow parameters at the inlet and outlet of the propulsor of C1 are the same as those of C2.


### 3.2. Boundary Conditions at Propulsor Inlet and Outlet

The propulsion inlet and outlet conditions are obtained using a component-based thermodynamic cycle analysis model of the propulsion system. It consists of a discrete thermodynamic cycle analysis module and an iterator. The thermodynamic cycle analysis module has three components, namely, diffuser, fan, and nozzle. The discrete Newton method is adopted as the iterator. The energy generation and transmission systems are not modeled in this model. The calculations are largely based on the formulation in the N+3 aircraft concept study report [8], and some modifications are made due to the different constitutions of the propulsion system. In this paper, we focus on the influence of BLI on the aerodynamic dissipation of the aircraft external flow; the efficiency variation of the propulsor with the flow distortion of BLI is out of our consideration. The decrease in stagnation pressure of the propulsor inflow caused by BLI is considered in the calculation, since the mass flow rate of the propulsor is highly affected by this parameter.The station numbers of the discrete thermodynamic cycle analysis module are shown in Figure6. Station 0 represents the free-stream condition, which is determined on the basis of the flight height and Mach number. Station 1 is the inlet plane where the fluid begins to enter the propulsor. Station 1.8 accounts for the stagnation pressure decrease due to boundary layer ingestion. Station 2, which is located at the fan face, is the starting point of pressure increase and the ending point of inlet flow. Station 2.1 is located behind the fan blades; it is the end point of the pressure increase and the starting point of the nozzle flow. Station 7 is the nozzle exit plane wherein fluid begins to leave the propulsor. Station 8 is the jet plane with ambient static pressure p∞. It is located a short distance behind the nozzle exit plane.Figure 6
Station numbers of the discrete thermodynamic cycle analysis module.The detailed constitution of the component-based propulsor model is shown in Figure7. From stations 0 to 1, the flow experienced a progression in isentropic compression, and 1D isentropic relations are used to obtain the flow parameters at station 1. A total pressure decrease parameter ΔPt is added between stations 1 and 1.8 to consider the kinematic energy loss of the boundary layer before the inlet. The flow path from stations 1.8 to 2 is the diffuser before the fan, which can be modeled by the total pressure recovery parameter πd. The flow path from stations 2 to 2.1 is the fan and modeled by total pressure ratio πf and polytropic efficiency ηf. The flow path from stations 2.1 to 7 is the nozzle, which is simulated by the total pressure recovery parameter πfn. An isentropic expansion is assumed from stations 7 to 8. The flow parameters at the jet plane (station 8) are obtained with the ambient static pressure condition (p=p∞). The discrete Newton iterator is used to update the fan mass flow rate mf and fan total pressure ratio πf. A fan operation map is utilized to obtain the polytropic efficiency ηf with given mf and πf. The fan operation map is generated from CFD calculations of the fan blades under different fan operation conditions. The constraint condition of the component-based propulsor model is a specified power requirement and a mass continuity constraint between the fan and nozzle. This study focuses on the influence of BLI on aircraft aerodynamic performance, and the influence of BLI on fan efficiency is ignored.Figure 7
Component-based propulsor model.The power requirement constraint Equation (26) is derived from the power balance Equation (8) that matches the propulsion system power supply with the power consumption under given flight conditions. That is, in this constraint, the propulsor should provide the net mechanical flow power to compensate for the specified airframe dissipation performance and the net streamwise force requirement.In addition to the mechanical power flow requirement, another constraint is that the mass flow rate of the fan and the nozzle should match. The constraint on mass continuity can be expressed as(37)mfρ∞V∞Sref−ρ7V7A7ρ∞V∞Sref=0,where the subscript7 indicates the properties in the nozzle exit plane. Quantity A7 denotes the exit plane area of the nozzle.


### 3.3. Postprogress Method

Calculating aerodynamic dissipation components from the CFD results requires the definition of the integral volume. As shown in Figure8, a cylinder containing the aircraft is constructed as the integral volume, and its axis is parallel to the free-stream velocity vector V∞. The nose of the aircraft is selected as the center of the integral cylinder. The radius of the cylinder rc and the distance from the aircraft trailing edge to the downstream boundary of the integral volume xd are controlled to show how the calculated dissipation components vary with the integration range. The upstream boundary is placed 10-folds of the root chord length croot, far from the aircraft nose where the flow parameters are equivalent to the free-stream parameters.Figure 8
Integral volume for aerodynamic dissipation calculations.The flow feature extraction method is required to attribute aerodynamic dissipation to various character flow regions and help establish a physical insightful understanding of the integrated airframe and propulsor performance. Fluid sensors are applied to the isolation of the jet, surface boundary layer, and shear layer regions. The trailing edge vortex dissipationΦvortex is calculated on the downstream boundary of the integral volume, according to Equation (17). Thus, a fluid sensor for extracting this flow feature is unnecessary.The sensor used to extract the propulsion jet [23] is
(38)γ=H−H∞H∞,whereH represents the total enthalpy of the isolated fluid and H∞ is the total enthalpy of free-stream fluid. The threshold of the jet scalar γ is selected as 0.04.The sensor for the shock wave extraction [23] is
(39)F=V⋅∇pa∇p,wherea is the local sound speed. The threshold of F is set to 0.95.The sensor for isolating the boundary layer and shear layer flow [23] is
(40)ξ=μ+μtμτ¯¯⋅∇⋅u⟶p∇⋅u⟶,whereμ is the molecular viscosity and μt is the turbulent viscosity. The threshold of shear layer scalar ξ is set to 1. This sensor is pragmatic given that viscous dissipation Φ is normalized by the absolutely volumetric pressure work PV to exclude the bulk strain rate contributions usually attributed to shocks.Integrating dissipations only in these extracted flow regions can omit the spurious contributions of the surrounding potential flow inΦ.


## 3.1. Aircraft Configuration Description

We consider a subsonic BWB transport aircraft designed by Yan et al. [27], which has a 6,000 nm air range with a cruise Mach of 0.85 and a cruise altitude of 35,000 ft. It can carry 350 passengers with a 507,495 lb maximum takeoff gross weight. The designed lift coefficient in the cruise condition is CL=0.36. Three different configurations are constructed, namely, C0, C1, and C2, to analyze the effects of aeropropulsor coupling on the aerodynamic performance of the BWB aircraft. C0 is a clean BWB airframe with no propulsor installed on it, which is actually the unpowered configuration. The BWB airframe is composed of three parts: center body, blend wing, and outer wing. C1 consists of a BWB airframe equivalent to C0 and a propulsor placed over the rear part of the upper surface of the center body of the airframe near the trailing edge. C1 represents the power-on configuration with a non-BLI propulsor. C2 has the same BWB airframe as C0 and C1 and a BLI propulsor embedded in the rear part of the upper surface of the center body of the airframe to ingest the boundary layer flow generated on the airframe surface. The detailed vertical and symmetry plane views of the three configurations are shown in Figure 5.Figure 5
Configurations used in the energy-based aerodynamic performance analysis.3D RANS calculations are performed for the three configurations. The free-stream Mach number is set to 0.85; and the Reynolds number is approximately7×107, which is the same under the aircraft cruise condition. Given that the BLI propulsion system is embedded in the rear part of the upper surface of the center body, the lift coefficients CL of C1 and C2 may differ from that of C0 at the same attack angle. Thus, we perform the calculation on a range of attack angles and compare the dissipation components of the three configurations with a similar lift coefficient and not the same angle. The initial values of flow parameters at the propulsor inlet and outlet of C2 are calculated using the component-based propulsor model described in Section 3.2. The aircraft dissipation components used in the component-based propulsor model are the results of the unpowered configuration C0. Then, a calculation matrix constructed by the attack angle array and the fan operation array is generated, which is used to find the cruise point (CL=0.36,CX=0) of C2. The flow parameters at the inlet and outlet of the propulsor of C1 are the same as those of C2.


## 3.2. Boundary Conditions at Propulsor Inlet and Outlet

The propulsion inlet and outlet conditions are obtained using a component-based thermodynamic cycle analysis model of the propulsion system. It consists of a discrete thermodynamic cycle analysis module and an iterator. The thermodynamic cycle analysis module has three components, namely, diffuser, fan, and nozzle. The discrete Newton method is adopted as the iterator. The energy generation and transmission systems are not modeled in this model. The calculations are largely based on the formulation in the N+3 aircraft concept study report [8], and some modifications are made due to the different constitutions of the propulsion system. In this paper, we focus on the influence of BLI on the aerodynamic dissipation of the aircraft external flow; the efficiency variation of the propulsor with the flow distortion of BLI is out of our consideration. The decrease in stagnation pressure of the propulsor inflow caused by BLI is considered in the calculation, since the mass flow rate of the propulsor is highly affected by this parameter.The station numbers of the discrete thermodynamic cycle analysis module are shown in Figure6. Station 0 represents the free-stream condition, which is determined on the basis of the flight height and Mach number. Station 1 is the inlet plane where the fluid begins to enter the propulsor. Station 1.8 accounts for the stagnation pressure decrease due to boundary layer ingestion. Station 2, which is located at the fan face, is the starting point of pressure increase and the ending point of inlet flow. Station 2.1 is located behind the fan blades; it is the end point of the pressure increase and the starting point of the nozzle flow. Station 7 is the nozzle exit plane wherein fluid begins to leave the propulsor. Station 8 is the jet plane with ambient static pressure p∞. It is located a short distance behind the nozzle exit plane.Figure 6
Station numbers of the discrete thermodynamic cycle analysis module.The detailed constitution of the component-based propulsor model is shown in Figure7. From stations 0 to 1, the flow experienced a progression in isentropic compression, and 1D isentropic relations are used to obtain the flow parameters at station 1. A total pressure decrease parameter ΔPt is added between stations 1 and 1.8 to consider the kinematic energy loss of the boundary layer before the inlet. The flow path from stations 1.8 to 2 is the diffuser before the fan, which can be modeled by the total pressure recovery parameter πd. The flow path from stations 2 to 2.1 is the fan and modeled by total pressure ratio πf and polytropic efficiency ηf. The flow path from stations 2.1 to 7 is the nozzle, which is simulated by the total pressure recovery parameter πfn. An isentropic expansion is assumed from stations 7 to 8. The flow parameters at the jet plane (station 8) are obtained with the ambient static pressure condition (p=p∞). The discrete Newton iterator is used to update the fan mass flow rate mf and fan total pressure ratio πf. A fan operation map is utilized to obtain the polytropic efficiency ηf with given mf and πf. The fan operation map is generated from CFD calculations of the fan blades under different fan operation conditions. The constraint condition of the component-based propulsor model is a specified power requirement and a mass continuity constraint between the fan and nozzle. This study focuses on the influence of BLI on aircraft aerodynamic performance, and the influence of BLI on fan efficiency is ignored.Figure 7
Component-based propulsor model.The power requirement constraint Equation (26) is derived from the power balance Equation (8) that matches the propulsion system power supply with the power consumption under given flight conditions. That is, in this constraint, the propulsor should provide the net mechanical flow power to compensate for the specified airframe dissipation performance and the net streamwise force requirement.In addition to the mechanical power flow requirement, another constraint is that the mass flow rate of the fan and the nozzle should match. The constraint on mass continuity can be expressed as(37)mfρ∞V∞Sref−ρ7V7A7ρ∞V∞Sref=0,where the subscript7 indicates the properties in the nozzle exit plane. Quantity A7 denotes the exit plane area of the nozzle.


## 3.3. Postprogress Method

Calculating aerodynamic dissipation components from the CFD results requires the definition of the integral volume. As shown in Figure8, a cylinder containing the aircraft is constructed as the integral volume, and its axis is parallel to the free-stream velocity vector V∞. The nose of the aircraft is selected as the center of the integral cylinder. The radius of the cylinder rc and the distance from the aircraft trailing edge to the downstream boundary of the integral volume xd are controlled to show how the calculated dissipation components vary with the integration range. The upstream boundary is placed 10-folds of the root chord length croot, far from the aircraft nose where the flow parameters are equivalent to the free-stream parameters.Figure 8
Integral volume for aerodynamic dissipation calculations.The flow feature extraction method is required to attribute aerodynamic dissipation to various character flow regions and help establish a physical insightful understanding of the integrated airframe and propulsor performance. Fluid sensors are applied to the isolation of the jet, surface boundary layer, and shear layer regions. The trailing edge vortex dissipationΦvortex is calculated on the downstream boundary of the integral volume, according to Equation (17). Thus, a fluid sensor for extracting this flow feature is unnecessary.The sensor used to extract the propulsion jet [23] is
(38)γ=H−H∞H∞,whereH represents the total enthalpy of the isolated fluid and H∞ is the total enthalpy of free-stream fluid. The threshold of the jet scalar γ is selected as 0.04.The sensor for the shock wave extraction [23] is
(39)F=V⋅∇pa∇p,wherea is the local sound speed. The threshold of F is set to 0.95.The sensor for isolating the boundary layer and shear layer flow [23] is
(40)ξ=μ+μtμτ¯¯⋅∇⋅u⟶p∇⋅u⟶,whereμ is the molecular viscosity and μt is the turbulent viscosity. The threshold of shear layer scalar ξ is set to 1. This sensor is pragmatic given that viscous dissipation Φ is normalized by the absolutely volumetric pressure work PV to exclude the bulk strain rate contributions usually attributed to shocks.Integrating dissipations only in these extracted flow regions can omit the spurious contributions of the surrounding potential flow inΦ.


## 4. Aerodynamic Performance Analysis of a BWB Aircraft with a BLI Propulsor

Firstly, the grid refinement, domain size, and integral volume size study of the calculated dissipation and force coefficients on the unpowered airframe configuration is displayed. Then, a quantitative evaluation of the BLI benefit applied on a subsonic BWB aircraft [27] and the D8 [21] aircraft is performed using the determining expression of PSC developed in Section 2.4. This evaluation of the benefits of the BLI is based on the dissipation components of the unpowered aircraft. Finally, the aerodynamic dissipation performance of the subsonic BWB aircraft with different propulsor installation types is analyzed. Through the comparison of the dissipation components of the BLI and non-BLI configurations, the aeropropulsive coupling effects can be expressed by the variation in the dissipation components and finally attributed to related flow features in the control volume, which provide a more physical insightful understanding of airframe–propulsor integration.The grid refinement, domain size, and integral volume size study is displayed in Section4.1. The quantitative evaluation of the BLI benefit is presented in Section 4.2. The detailed analysis of the aeropropulsive coupling effects of the BWB aircraft with a BLI propulsor is presented in Section 4.3.


### 4.1. Grid Refinement, Domain Size, and Integral Volume Size Study of Unpowered ConfigurationC0

The clean-airframe configurationC0 is chosen to conduct the grid refinement, domain size, and integral volume size study. The front two part is used to exhibit the variation in the lift, drag, and dissipation coefficients with the cell number and domain size of the mesh. Then, the integral volume size study is used to find a proper size of the integral control volume and a proper position of the TP plane for the calculation of dissipation components. As shown in Figure 9, the meshes constructed here are C-type meshes, and the distance of downstream wake plane ld from the airframe trailing edge is three times that of other outer boundaries lΣ from the airframe surface. The definition of the integral control volume used for dissipation component calculation is previously presented in Figure 8.Figure 9
Definition of mesh domain used in CFD calculations.The liftCL and drag CD coefficients are referenced to 1/2ρ∞V∞2Sref. The total aerodynamic dissipation coefficient CΦ∗ is referenced to 1/2ρV∞3Sref, where Sref is the reference area. Calculations are performed at Reynolds number Re∞=7×107 and free-stream Mach number M∞=0.85. The angle of attack is AOA=3.8°. CL and CD are obtained by integrating the pressure and viscous forces on the configuration surface of C0 and then projecting the force vector to the specialized direction. CΦ∗ is obtained using Equation (10). In the integral volume study, CD serves as the reference value for CΦ∗. The difference between CΦ∗ and CD is an indicator of the accuracy of the calculation of CΦ∗.


#### 4.1.1. Grid Refinement Study

Three meshes with different cell numbers are constructed aroundC0 to carry out the grid refinement study. The outer boundary of the meshes is 30croot, far from the airframe surface. The height of the first cell layer from the surface of the aircraft is set to 5×10−6m to ensure y+<1. The following CFD results show that y+≈0.4. The cell numbers of the three meshes are 3.80 million, 7.73 million, and 10.73 million. The detailed setup of the grid refinement study is shown in Table 1.Table 1
Grid refinement study setup.
Cell numberhy+lΣ/crootCoarse3.80 million3.100.430crootMedium7.73 million1.520.430crootFine10.73 million1.070.430crootAn indicator of the size of the mesh cellh is defined as
(41)h=1N∑i=1NΔVi1/3,whereN is the cell number involved in the integral calculation of CΦ∗.Table2 shows that the calculated CL, CD, and CΦ∗ have a good convergence characteristic with the mesh refinement. The maximum variation in the calculated CL is below 1% of the value of the medium mesh with h=1.07. The integral results of CD show that the difference between the values of the coarse and the medium mesh is approximately 7% of the value of the medium mesh, while the difference between the values of the medium and the fine mesh is merely 0.6% of the value of the fine mesh. The calculated CΦ∗ results show that the difference in the values of the coarse and medium meshes in CΦ∗ is nearly 1% of the value of the medium mesh, and the CΦ∗ difference between the medium and the fine mesh is also 1% of the value of the fine mesh. The maximum difference in CD and CΦ∗ appears on the coarse mesh, which is around 9% of the value of CD. The minimum difference in CD and CΦ∗ appears on the fine mesh, which is approximately 2% of the value of CD. The results of CΦ∗ is insensitive to the mesh cell number if y+<1 is satisfied and the boundary layer and wake regions are contained in the integral volume, which shows the advantage of the energy-based dissipation coefficient CΦ∗ over the momentum-based drag coefficient CD in the aerodynamic performance evaluation of unpowered aircraft.Table 2
Grid refinement study result.
hCLCDCΦ∗Coarse3.100.3710.01950.0177Medium1.520.3730.01820.0175Fine1.070.3740.01810.0177


#### 4.1.2. Domain Size Study

Three meshes with different domain sizes are constructed aroundC0 to conduct the domain size study. The outer boundaries of the meshes are 20croot, 30croot, and 40croot, far from the surface of the airframe, where croot is the root airfoil chord length of the airframe. The height of the first cell layer from the surface of the aircraft is set to 5×10−6m to ensure y+<1. The following CFD results show that y+≈0.4. The cell numbers of the three meshes are 6.99 million, 7.73 million, and 8.04 million. The domain size is represented by a nondimensionalized indicator lΣ/croot. The detailed setup of the domain size study is shown in Table 3.Table 3
Domain size study setup.
Cell numbery+lΣ/crootD16.99 million0.420crootD27.73 million0.430crootD38.04 million0.440crootTable4 shows that the calculated CL, CD, and CΦ∗ are insensitive to the mesh domain size. The maximum variation in CL is approximately 0.3% of the value of the D2 mesh (the same as the medium mesh in the mesh refinement study) withlΣ/croot. The maximum variation in CD is nearly 0.5% of the value of the D2 mesh. The difference in CΦ∗ between D1 and D2 meshes is about 1% of the value of the D2 mesh. The difference in CΦ∗ between D2 and D3 meshes is similar to that between D1 and D2. The differences in CΦ∗ and CD of the same mesh are 3%, 4%, and 5% for D3, D2, and D1 meshes, respectively. CΦ∗ is more sensitive to the size of the mesh domain than CD and CL.Table 4
Domain size study result.
lΣ/crootCLCDCΦ∗D120croot0.3730.01820.0173D230croot0.3730.01820.0175D340croot0.3740.01830.0177


#### 4.1.3. Integral Volume Size Study

The artificial dissipation introduced by the numerical simulation brings spurious effects on the calculation of dissipation components. Thus, the choice of the integral domain size is important for the prediction of the dissipation components. Usually, the drag coefficientCD calculated on the aircraft surface is chosen as the reference value to evaluate the accuracy of the total dissipation coefficient calculation. The variations in the calculated total dissipation and dissipation components of C0 with the integral volume size indicator rc and xd are analyzed. The results are presented in Figure 10. All the calculated dissipation coefficients are normalized by nearfield drag coefficient CD. The integral volume size indicators rc and xd are normalized by croot.Figure 10
Sensitivity ofCΦ∗, CΦ, −CPV, and CĖ to the integral volume size. The values are normalized by near-field drag coefficient CD. (a) Variation in CΦ∗, CΦ, −CPV, and CĖ with normalized integral volume size indicator rc/croot. (b) Variation in CΦ∗, CΦ, −CPV, and CĖ with normalized integral volume size indicator xd/croot.
(a)(b)The variations inCΦ∗, CΦ, −CPV, and CĖ with rc/croot are presented in Figure 10(a), with a fixed xd/croot≈1.24. The total dissipation coefficients CΦ∗ increase with the integral domain size indicator rc/croot, which is a cumulative result of the increase in CΦ and CĖ with rc/croot. The dissipation component −CPV does not change with rc/croot when rc/croot>1 because the static pressure is recovered to ambient pressure at this location. The proper range of rc/croot is 2–4, in which the difference between calculated CΦ∗ and the near-field drag coefficient CD is within 5% of the value of CD.The variations inCΦ∗, CΦ, −CPV, and CĖ with xd/croot are presented in Figure 10(b), with a fixed rc/croot≈2.24. The calculated total dissipation CΦ∗ increases with xd/croot because of numerical spurious effects. When moving the TP (downstream outer boundary of the integral volume) downstream, energy transfer occurs among CΦ, −CPV, and CĖ. The volumetric pressure work coefficient −CPV represents the rate of energy transferred from mechanical power to internal energy, most of which are reversible. The decrease in −CPV with xd/croot means that the inner energy of the fluid is transferred to the mechanical energy of the fluid through an expansion progress, which is completed at xd/croot=1.24. When xd/croot>1.24, coefficient −CPV does not change with xd/croot anymore. xd/croot=1.24 is a proper place for the placement of TP, where the reversible energy transferring progress is completed and the calculated −CPV and CĖ do not contain the interference of the potential flow region. The remaining −CPV is about 13% of the total dissipation coefficient CΦ∗, which may indicate there are shock wave regions in the flow field. This term will be discussed next.Energy transfer occurs among decomposed components ofCĖ. The variation in the decomposed constituents of CĖ with xd/croot is shown in Figure 11. At the trailing edge of the aircraft (xd/croot=0), the absolute values of CĖa, CĖv, and CĖp are all over 50% of CD, while the overall CĖ sum is only 10% of CD. The signs of CĖa and CĖv are contrary to that of CĖp. Thus, the large absolute values of CĖa, CĖv, and CĖp are canceled in CĖ. The large absolute values of CĖa, CĖv, and CĖp near the aircraft trailing edge are caused by the potential flow region interference, which quickly decays within a fraction of croot.Figure 11
Sensitivity ofCĖa, CĖv, CĖp, and CĖ to the streamwise position of TP. The values of calculated dissipation components are normalized by near-field drag coefficient CD.Figures10 and 11 show the progress in energy transfer between the decomposed components of CΦ∗ and CĖ. Based on the variation charterer of −CPV and decomposed constitutions of CĖ, xd/croot=1.24 is a proper place for the placement of TP for the decomposition of CΦ∗, where the reversible energy transferring progress is completed and the calculated −CPV and CĖ components do not contain the interference of the potential flow region. At xd/croot=1.24, CĖa, and CĖp are nearly decaying to zero, and CĖ only accounts for the energy deposited in the trailing edge vortex, which is defined as the vortex dissipation component. All the viscous dissipation that occurred before xd/croot=1.24 was defined as the profile dissipation component. The no longer changing −CPV is treated as the volume pressure work-related dissipation components, which will be nonzero when the irreversible compression–expansion progress occurs [23].


### 4.2. Evaluation of BLI Benefits and Its Limitation of the BWB Aircraft with a BLI Propulsor

The power balance analysis framework presented in Section2.4. is used to quantify the BLI benefit of the BWB aircraft [27] and the D8 aircraft [21]. The results of the dissipation decomposition of the unpowered configuration are used in the analysis of the BLI benefits. The values of parameters fBLI, α, β, and ε, defined by Equations (24), (25), (27), and (29), are shown in Table 5, together with those of PSC, PSCmax, PSCjet, and PSCaf.Table 5
Comparison of PSC for the BWB aircraft and D8 [21].
fBLIαβεPSCPSCmaxPSCjetPSCafBWB22.5%9.0%60%0.35%18.8%3%2%D817.0%8.5%65%0.66%26.8%5%1%As shown in Table5, The power-saving coefficient PSC of the BWB aircraft [27] and D8 [21] is 5% and 6%, respectively. The evaluated limitations of the BLI benefit PSCmax of the two BLI aircrafts are 18.8% and 26.8%, respectively; this range is similar to the experiment results of Lv et al. [19, 20]. Despite having a higher value of fBLI, the value of PSC in the BWB aircraft [27] is still lower than that in D8 [21], which is mainly due to a lower value of ε.According to Equations (35) and (36), the term for the airframe dissipation reduction PSCaf and the term for the jet dissipation reduction PSCjet are largely determined by the propulsor design parameter ε. The value of PSCaf decreases with ε, while the value of PSCjet increases with ε. Thus, given the similar values of α and β, PSCjet is more dominant and PSCaf is weaker when the value of ε is higher. As indicated in Table 5, the value of ε for D8 is one-time greater than that for the BWB aircraft. As a result, the values of PSCjet and PSCaf are nearly the same on the BWB aircraft while in the PSC constitution of the D8 aircraft, PSCjet is the dominant factor. Additionally, the PSCjet of D8 is more evident than that of the BWB aircraft; that is, it is nearly one-time higher.


### 4.3. Energy-Based Aerodynamic Analysis of Aeropropulsion Coupling Effects

TheC2 case near the cruise point (CL=0.36, CX=0) is selected to conduct the dissipation decomposition. The cases of C0 and C1 with nearly a similar lift coefficient are chosen to conduct the dissipation constitution comparison with C2. In the CFD postprogress, PK is obtained from the mechanical flow power integral in Equation (3) on the inlet and outlet planes of the propulsion system, FX is obtained by integrating pressure and viscous force axial projection onto the airframe’s surface and axial momentum flux across the inlet and outlet planes, and lift force L is calculated by integrating pressure and viscous force vertical axial projection over the airframe’s surface. Total aerodynamic dissipation Φ∗ can be indirectly obtained through the power balance equation, Equation (8), or be calculated directly from the viscous dissipation Φ, volumetric pressure work PV, and the irreversible mechanical energy outflow Ė integral in Equation (7). The difference between the calculated Φ∗ from Equations (7) and (8) is within 5% of the value calculated from Equation (8) for all the three configurations. The vortex dissipation Φvortex in Equations (14) and ((17)) are calculated by the integral of Ė on the TP, according to Equation (17). The lift CL and net streamwise force CX coefficients are referenced to 1/2ρ∞V∞2Sref. The energy-related coefficients CPK, CΦ∗, CĖ, CPV, and CΦ are referenced to 1/2ρV∞3Sref, where Sref is the reference area.The results of the dissipation decomposition using Equation (7) are presented in Table 6. The ratio of CΦ to CΦ∗ of the three configurations is about 60%~70%, which contains the contribution from the aircraft surface boundary layer, trailing edge wake, and jet. The ratio of −CPV to CΦ∗ is quite considerable (10%~15%), which indicates there are shock wave regions in the flow field. The ratio of CĖ to CΦ∗ is about 18%~26%, lower than that of CΦ to CΦ∗.Table 6
Dissipation components forC0, C1, and C2.
CΦCĖ−CPVCΦ∗C00.0107(60.5%)0.0046(26.0%)0.0024(13.5%)0.0177C10.0190(66.7%)0.0052(18.2%)0.0043(15.1%)0.0285C20.0138(67.0%)0.0049(23.8%)0.(9.2%)0.0206To explain the variation inCΦ, the flow feature extraction methods described in Section 3.3. are used to isolate the jet region, the shock wave region, and the boundary layer/wake region in the flow field. The extraction priority is first the jet region, then the shock wave region, and last the shear layer region. The extracted regions are shown in Figure 12. The calculated viscous dissipation coefficients in the extracted regions are shown in Figure 13.Figure 12
Jet region, shock region, and shear layer region.Figure 13
Dissipation component comparison ofC0, C1, and C2.As is shown in Figure13, it is the difference in the dissipation that occurred in the shear layer region that causes the difference in CΦ. Then, the dissipation coefficient of shear layer CΦp is divided into five parts using the position-related sensor to show the variation in the dissipation generated on different parts of the BWB aircraft surfaces and in the wakes. The positions of the five dissipation components are labeled in Figure 14. Φcb‐bl represents the dissipation that occurred on the surfaces of the center body. Φow represents the dissipation that occurred on the surfaces of the blend wing and outer wing. Φeng‐bl represents the dissipation that occurred on the surfaces of the propulsor nacelle. Φcb‐wake represents the dissipation that occurred in the center body trailing edge wake. Φeng‐wake represents the dissipation that occurred in the propulsor nacelle trailing edge wake.Figure 14
The position ofΦcb‐bl, Φow, Φeng‐bl, Φcb‐wake, and Φeng‐wake.
(a)(b)According to the results of the decomposition of the shear layer dissipation coefficientCΦp of C0, C1, and C2 in Table 7, the dissipation components with obvious variation are CΦcb‐wake, CΦeng‐bl, and CΦeng‐wake. By comparing the CΦcb‐wake term of C1 and C2, the wake dissipation reduction brought by BLI is 0.0005, while that calculated from Equations (36) and (32) is 0.0004. These results show that the airframe dissipation reduction is underestimated by Equation (36), mainly due to the fact that we ignore the change in the velocity profile of the trailing edge of the center body due to propulsor installation. Meanwhile, the variations in CΦcb‐bl and CΦow of the three configurations are less obvious. The value of CΦcb‐bl of C2 is slightly less than that of C0 and C1 due to the reduction in the surface area of the center body. The large increase in CΦeng‐bl and CΦeng‐wake may be due to the flow separation on the propulsor outer surface as is shown in Figure 15, which can be reduced through a parametric design of the configuration.Table 7
Shear layer dissipation coefficientCΦp decomposition.
CΦcb‐blCΦcb‐wakeCΦeng‐blCΦeng‐wakeCΦowCΦpC00.00520.0004000.00450.0101C10.00520.00100.00330.00160.00440.0155C20.00480.00050.00080.00020.00440.0107Figure 15
The contours of Mach number and streamline ofC1.To explain the variation inCΦcb‐wake and CΦeng‐wake, the trailing edge velocity profiles of the three configurations (C0, C1, and C2) have been extracted from the previously obtained CFD results, which are presented in Figure 16. The extraction positions of the trailing edge velocity profiles are marked in Figures 16(a)–16(c). Figures 16(d)–16(f) are the trailing edge velocity profiles of the center body and propulsor for C0, C1, and C2, respectively. The center body trailing edge velocity profile of C1 (Figure 16(e)) is steeper than that of C0 (Figure 16(d)), which represents a larger velocity gradient of the boundary layer flow at the trailing edge and thus a larger CΦcb‐wake. The center body trailing edge velocity profile of C2 (Figure 16(f)) is moderate compared with that of C1 (Figure 16(e)), which is the result of the wake filling effect of BLI. The moderate center body trailing edge velocity profile of C2 represents a smaller velocity gradient of the boundary layer flow and brings a reduction in CΦcb‐wake. The reduction in CΦeng‐wake of C2 compared with C1 can also be attributed to a moderated trailing edge velocity profile.Figure 16
Trailing edge velocity profiles forC0, C1, and C2. (a–c) Positions of the extracted trailing edge velocity profiles for C0, C1, and C2. (d–f) Trailing edge velocity profiles of center body and propulsor for C0, C1, and C2.
(a)(b)(c)(d)(e)(f)The volumetric pressure work coefficient−CPV of C0, C1, and C2 is 0.0024, 0.0041, and 0.0019, respectively. The value of −CPV of C2 is smaller than that of C1, which means a reduction in the waste of mechanical energy. As is explained by Sanders and Laskaridis [23], the value of −CPV represents the mechanical energy waste in the control volume that results in an increase in internal energy. This term cannot be neglected when there is a shock wave in the flow field. When the fluid passes through a shock wave region, the mechanical energy will decrease because of expansion damping [23].In order to show the influence of the shock wave on the variation of−CPV, the flow fields at Ma=0.85 and Ma=0.6 of C0 were generated by the CFD method. Then, the 2D slices of the two flow fields at the spanwise location Z=2m were extracted, and the contours of velocity divergence ∇⋅V and relative pressure contours p−p∞ are shown in Figure 17. In this figure, the contour flood is for ∇⋅V, and the contour line is for p−p∞. Figure 17(a) is extracted at Ma=0.85, and Figure 17(b) is extracted at Ma=0.6. As is shown in Figure 17(a), there is a strong compression region in the high Mach number (Ma=0.85) flow, identified by the large negative value of ∇⋅V. This strong compression region is actually the shock wave region, and there may be obvious mechanical energy loss in this region. Thus, the contribution of −CPV to CΦ∗ cannot be neglected. However, in Figure 17(b), there is no obvious compression region in the flow field with Ma=0.6. Thus, the contribution of −CPV to CΦ∗ is only 4%, which can be neglected. The above comparison shows that −CPV is highly affected by the existence of the shock wave region and, thus, is highly related to the mechanical energy loss due to the shock wave. The obvious variation in −CPV among C0, C1, and C2 may due to the change of strength of the shock wave.Figure 17
The contour of velocity divergence∇⋅V and relative pressure p−p∞. The contour flood is for ∇⋅V, and the contour line is for p−p∞. The labeled figures are the value of p−p∞. (a) The contours at Ma=0.85. (b) The contours at Ma=0.6.
(a)(b)The coefficients of the dissipations that occurred in the trailing edge vortexCĖ of C0, C1, and C2 are 0.0046, 0.0052, and 0.0049, respectively. The variation in CĖ between C1 and C2 shows that the installation height of the propulsor system can change the lift coefficient distribution along the aircraft spanwise, which is the determining term of lift-induced vortex strength. To explain the variation in CĖ of C0, C1, and C2, the contours of Ev_uni on TP located at X=100m (2.3croot) of C0, C1, and C2 are generated from the CFD results, as shown in Figure 18. Ev_uni is the nondimensionalized transverse kinetic energy deposited in the trailing edge vortex, and it is defined as
(42)Ev_uni=1/2v2+w21/2V∞2.Figure 18
Ev_uni contours of C0, C1, and C2. (a–c) TP plane positions. (d–f) Ev_uni contours on TP plane.
(a)(b)(c)(d)(e)(f)The dark red region in Figure18(e) represents a notable contribution to Φvortex of C1, which is weaker in the Ev_uni contour of C0 and C2 (Figures 18(d)–18(f)). The difference in Ev_uni distribution of C0, C1, and C2 shows the source of variation in Φvortex in all three configurations.Summarizing the abovementioned analysis on the results of the total dissipation coefficientCΦ∗ breakdown of C0, C1, and C2 for the BWB aircraft shows that the airframe dissipation components that are strongly affected by BLI effect are −CPV, CΦwake, and CΦvortex. Variations in CΦeng‐bl and CΦeng‐wake are not strictly caused by BLI effects, but they are a result of inappropriate design of the propulsor nacelle that can be improved later.


## 4.1. Grid Refinement, Domain Size, and Integral Volume Size Study of Unpowered ConfigurationC0

The clean-airframe configurationC0 is chosen to conduct the grid refinement, domain size, and integral volume size study. The front two part is used to exhibit the variation in the lift, drag, and dissipation coefficients with the cell number and domain size of the mesh. Then, the integral volume size study is used to find a proper size of the integral control volume and a proper position of the TP plane for the calculation of dissipation components. As shown in Figure 9, the meshes constructed here are C-type meshes, and the distance of downstream wake plane ld from the airframe trailing edge is three times that of other outer boundaries lΣ from the airframe surface. The definition of the integral control volume used for dissipation component calculation is previously presented in Figure 8.Figure 9
Definition of mesh domain used in CFD calculations.The liftCL and drag CD coefficients are referenced to 1/2ρ∞V∞2Sref. The total aerodynamic dissipation coefficient CΦ∗ is referenced to 1/2ρV∞3Sref, where Sref is the reference area. Calculations are performed at Reynolds number Re∞=7×107 and free-stream Mach number M∞=0.85. The angle of attack is AOA=3.8°. CL and CD are obtained by integrating the pressure and viscous forces on the configuration surface of C0 and then projecting the force vector to the specialized direction. CΦ∗ is obtained using Equation (10). In the integral volume study, CD serves as the reference value for CΦ∗. The difference between CΦ∗ and CD is an indicator of the accuracy of the calculation of CΦ∗.


### 4.1.1. Grid Refinement Study

Three meshes with different cell numbers are constructed aroundC0 to carry out the grid refinement study. The outer boundary of the meshes is 30croot, far from the airframe surface. The height of the first cell layer from the surface of the aircraft is set to 5×10−6m to ensure y+<1. The following CFD results show that y+≈0.4. The cell numbers of the three meshes are 3.80 million, 7.73 million, and 10.73 million. The detailed setup of the grid refinement study is shown in Table 1.Table 1
Grid refinement study setup.
Cell numberhy+lΣ/crootCoarse3.80 million3.100.430crootMedium7.73 million1.520.430crootFine10.73 million1.070.430crootAn indicator of the size of the mesh cellh is defined as
(41)h=1N∑i=1NΔVi1/3,whereN is the cell number involved in the integral calculation of CΦ∗.Table2 shows that the calculated CL, CD, and CΦ∗ have a good convergence characteristic with the mesh refinement. The maximum variation in the calculated CL is below 1% of the value of the medium mesh with h=1.07. The integral results of CD show that the difference between the values of the coarse and the medium mesh is approximately 7% of the value of the medium mesh, while the difference between the values of the medium and the fine mesh is merely 0.6% of the value of the fine mesh. The calculated CΦ∗ results show that the difference in the values of the coarse and medium meshes in CΦ∗ is nearly 1% of the value of the medium mesh, and the CΦ∗ difference between the medium and the fine mesh is also 1% of the value of the fine mesh. The maximum difference in CD and CΦ∗ appears on the coarse mesh, which is around 9% of the value of CD. The minimum difference in CD and CΦ∗ appears on the fine mesh, which is approximately 2% of the value of CD. The results of CΦ∗ is insensitive to the mesh cell number if y+<1 is satisfied and the boundary layer and wake regions are contained in the integral volume, which shows the advantage of the energy-based dissipation coefficient CΦ∗ over the momentum-based drag coefficient CD in the aerodynamic performance evaluation of unpowered aircraft.Table 2
Grid refinement study result.
hCLCDCΦ∗Coarse3.100.3710.01950.0177Medium1.520.3730.01820.0175Fine1.070.3740.01810.0177


### 4.1.2. Domain Size Study

Three meshes with different domain sizes are constructed aroundC0 to conduct the domain size study. The outer boundaries of the meshes are 20croot, 30croot, and 40croot, far from the surface of the airframe, where croot is the root airfoil chord length of the airframe. The height of the first cell layer from the surface of the aircraft is set to 5×10−6m to ensure y+<1. The following CFD results show that y+≈0.4. The cell numbers of the three meshes are 6.99 million, 7.73 million, and 8.04 million. The domain size is represented by a nondimensionalized indicator lΣ/croot. The detailed setup of the domain size study is shown in Table 3.Table 3
Domain size study setup.
Cell numbery+lΣ/crootD16.99 million0.420crootD27.73 million0.430crootD38.04 million0.440crootTable4 shows that the calculated CL, CD, and CΦ∗ are insensitive to the mesh domain size. The maximum variation in CL is approximately 0.3% of the value of the D2 mesh (the same as the medium mesh in the mesh refinement study) withlΣ/croot. The maximum variation in CD is nearly 0.5% of the value of the D2 mesh. The difference in CΦ∗ between D1 and D2 meshes is about 1% of the value of the D2 mesh. The difference in CΦ∗ between D2 and D3 meshes is similar to that between D1 and D2. The differences in CΦ∗ and CD of the same mesh are 3%, 4%, and 5% for D3, D2, and D1 meshes, respectively. CΦ∗ is more sensitive to the size of the mesh domain than CD and CL.Table 4
Domain size study result.
lΣ/crootCLCDCΦ∗D120croot0.3730.01820.0173D230croot0.3730.01820.0175D340croot0.3740.01830.0177


### 4.1.3. Integral Volume Size Study

The artificial dissipation introduced by the numerical simulation brings spurious effects on the calculation of dissipation components. Thus, the choice of the integral domain size is important for the prediction of the dissipation components. Usually, the drag coefficientCD calculated on the aircraft surface is chosen as the reference value to evaluate the accuracy of the total dissipation coefficient calculation. The variations in the calculated total dissipation and dissipation components of C0 with the integral volume size indicator rc and xd are analyzed. The results are presented in Figure 10. All the calculated dissipation coefficients are normalized by nearfield drag coefficient CD. The integral volume size indicators rc and xd are normalized by croot.Figure 10
Sensitivity ofCΦ∗, CΦ, −CPV, and CĖ to the integral volume size. The values are normalized by near-field drag coefficient CD. (a) Variation in CΦ∗, CΦ, −CPV, and CĖ with normalized integral volume size indicator rc/croot. (b) Variation in CΦ∗, CΦ, −CPV, and CĖ with normalized integral volume size indicator xd/croot.
(a)(b)The variations inCΦ∗, CΦ, −CPV, and CĖ with rc/croot are presented in Figure 10(a), with a fixed xd/croot≈1.24. The total dissipation coefficients CΦ∗ increase with the integral domain size indicator rc/croot, which is a cumulative result of the increase in CΦ and CĖ with rc/croot. The dissipation component −CPV does not change with rc/croot when rc/croot>1 because the static pressure is recovered to ambient pressure at this location. The proper range of rc/croot is 2–4, in which the difference between calculated CΦ∗ and the near-field drag coefficient CD is within 5% of the value of CD.The variations inCΦ∗, CΦ, −CPV, and CĖ with xd/croot are presented in Figure 10(b), with a fixed rc/croot≈2.24. The calculated total dissipation CΦ∗ increases with xd/croot because of numerical spurious effects. When moving the TP (downstream outer boundary of the integral volume) downstream, energy transfer occurs among CΦ, −CPV, and CĖ. The volumetric pressure work coefficient −CPV represents the rate of energy transferred from mechanical power to internal energy, most of which are reversible. The decrease in −CPV with xd/croot means that the inner energy of the fluid is transferred to the mechanical energy of the fluid through an expansion progress, which is completed at xd/croot=1.24. When xd/croot>1.24, coefficient −CPV does not change with xd/croot anymore. xd/croot=1.24 is a proper place for the placement of TP, where the reversible energy transferring progress is completed and the calculated −CPV and CĖ do not contain the interference of the potential flow region. The remaining −CPV is about 13% of the total dissipation coefficient CΦ∗, which may indicate there are shock wave regions in the flow field. This term will be discussed next.Energy transfer occurs among decomposed components ofCĖ. The variation in the decomposed constituents of CĖ with xd/croot is shown in Figure 11. At the trailing edge of the aircraft (xd/croot=0), the absolute values of CĖa, CĖv, and CĖp are all over 50% of CD, while the overall CĖ sum is only 10% of CD. The signs of CĖa and CĖv are contrary to that of CĖp. Thus, the large absolute values of CĖa, CĖv, and CĖp are canceled in CĖ. The large absolute values of CĖa, CĖv, and CĖp near the aircraft trailing edge are caused by the potential flow region interference, which quickly decays within a fraction of croot.Figure 11
Sensitivity ofCĖa, CĖv, CĖp, and CĖ to the streamwise position of TP. The values of calculated dissipation components are normalized by near-field drag coefficient CD.Figures10 and 11 show the progress in energy transfer between the decomposed components of CΦ∗ and CĖ. Based on the variation charterer of −CPV and decomposed constitutions of CĖ, xd/croot=1.24 is a proper place for the placement of TP for the decomposition of CΦ∗, where the reversible energy transferring progress is completed and the calculated −CPV and CĖ components do not contain the interference of the potential flow region. At xd/croot=1.24, CĖa, and CĖp are nearly decaying to zero, and CĖ only accounts for the energy deposited in the trailing edge vortex, which is defined as the vortex dissipation component. All the viscous dissipation that occurred before xd/croot=1.24 was defined as the profile dissipation component. The no longer changing −CPV is treated as the volume pressure work-related dissipation components, which will be nonzero when the irreversible compression–expansion progress occurs [23].


## 4.1.1. Grid Refinement Study

Three meshes with different cell numbers are constructed aroundC0 to carry out the grid refinement study. The outer boundary of the meshes is 30croot, far from the airframe surface. The height of the first cell layer from the surface of the aircraft is set to 5×10−6m to ensure y+<1. The following CFD results show that y+≈0.4. The cell numbers of the three meshes are 3.80 million, 7.73 million, and 10.73 million. The detailed setup of the grid refinement study is shown in Table 1.Table 1
Grid refinement study setup.
Cell numberhy+lΣ/crootCoarse3.80 million3.100.430crootMedium7.73 million1.520.430crootFine10.73 million1.070.430crootAn indicator of the size of the mesh cellh is defined as
(41)h=1N∑i=1NΔVi1/3,whereN is the cell number involved in the integral calculation of CΦ∗.Table2 shows that the calculated CL, CD, and CΦ∗ have a good convergence characteristic with the mesh refinement. The maximum variation in the calculated CL is below 1% of the value of the medium mesh with h=1.07. The integral results of CD show that the difference between the values of the coarse and the medium mesh is approximately 7% of the value of the medium mesh, while the difference between the values of the medium and the fine mesh is merely 0.6% of the value of the fine mesh. The calculated CΦ∗ results show that the difference in the values of the coarse and medium meshes in CΦ∗ is nearly 1% of the value of the medium mesh, and the CΦ∗ difference between the medium and the fine mesh is also 1% of the value of the fine mesh. The maximum difference in CD and CΦ∗ appears on the coarse mesh, which is around 9% of the value of CD. The minimum difference in CD and CΦ∗ appears on the fine mesh, which is approximately 2% of the value of CD. The results of CΦ∗ is insensitive to the mesh cell number if y+<1 is satisfied and the boundary layer and wake regions are contained in the integral volume, which shows the advantage of the energy-based dissipation coefficient CΦ∗ over the momentum-based drag coefficient CD in the aerodynamic performance evaluation of unpowered aircraft.Table 2
Grid refinement study result.
hCLCDCΦ∗Coarse3.100.3710.01950.0177Medium1.520.3730.01820.0175Fine1.070.3740.01810.0177


## 4.1.2. Domain Size Study

Three meshes with different domain sizes are constructed aroundC0 to conduct the domain size study. The outer boundaries of the meshes are 20croot, 30croot, and 40croot, far from the surface of the airframe, where croot is the root airfoil chord length of the airframe. The height of the first cell layer from the surface of the aircraft is set to 5×10−6m to ensure y+<1. The following CFD results show that y+≈0.4. The cell numbers of the three meshes are 6.99 million, 7.73 million, and 8.04 million. The domain size is represented by a nondimensionalized indicator lΣ/croot. The detailed setup of the domain size study is shown in Table 3.Table 3
Domain size study setup.
Cell numbery+lΣ/crootD16.99 million0.420crootD27.73 million0.430crootD38.04 million0.440crootTable4 shows that the calculated CL, CD, and CΦ∗ are insensitive to the mesh domain size. The maximum variation in CL is approximately 0.3% of the value of the D2 mesh (the same as the medium mesh in the mesh refinement study) withlΣ/croot. The maximum variation in CD is nearly 0.5% of the value of the D2 mesh. The difference in CΦ∗ between D1 and D2 meshes is about 1% of the value of the D2 mesh. The difference in CΦ∗ between D2 and D3 meshes is similar to that between D1 and D2. The differences in CΦ∗ and CD of the same mesh are 3%, 4%, and 5% for D3, D2, and D1 meshes, respectively. CΦ∗ is more sensitive to the size of the mesh domain than CD and CL.Table 4
Domain size study result.
lΣ/crootCLCDCΦ∗D120croot0.3730.01820.0173D230croot0.3730.01820.0175D340croot0.3740.01830.0177


## 4.1.3. Integral Volume Size Study

The artificial dissipation introduced by the numerical simulation brings spurious effects on the calculation of dissipation components. Thus, the choice of the integral domain size is important for the prediction of the dissipation components. Usually, the drag coefficientCD calculated on the aircraft surface is chosen as the reference value to evaluate the accuracy of the total dissipation coefficient calculation. The variations in the calculated total dissipation and dissipation components of C0 with the integral volume size indicator rc and xd are analyzed. The results are presented in Figure 10. All the calculated dissipation coefficients are normalized by nearfield drag coefficient CD. The integral volume size indicators rc and xd are normalized by croot.Figure 10
Sensitivity ofCΦ∗, CΦ, −CPV, and CĖ to the integral volume size. The values are normalized by near-field drag coefficient CD. (a) Variation in CΦ∗, CΦ, −CPV, and CĖ with normalized integral volume size indicator rc/croot. (b) Variation in CΦ∗, CΦ, −CPV, and CĖ with normalized integral volume size indicator xd/croot.
(a)(b)The variations inCΦ∗, CΦ, −CPV, and CĖ with rc/croot are presented in Figure 10(a), with a fixed xd/croot≈1.24. The total dissipation coefficients CΦ∗ increase with the integral domain size indicator rc/croot, which is a cumulative result of the increase in CΦ and CĖ with rc/croot. The dissipation component −CPV does not change with rc/croot when rc/croot>1 because the static pressure is recovered to ambient pressure at this location. The proper range of rc/croot is 2–4, in which the difference between calculated CΦ∗ and the near-field drag coefficient CD is within 5% of the value of CD.The variations inCΦ∗, CΦ, −CPV, and CĖ with xd/croot are presented in Figure 10(b), with a fixed rc/croot≈2.24. The calculated total dissipation CΦ∗ increases with xd/croot because of numerical spurious effects. When moving the TP (downstream outer boundary of the integral volume) downstream, energy transfer occurs among CΦ, −CPV, and CĖ. The volumetric pressure work coefficient −CPV represents the rate of energy transferred from mechanical power to internal energy, most of which are reversible. The decrease in −CPV with xd/croot means that the inner energy of the fluid is transferred to the mechanical energy of the fluid through an expansion progress, which is completed at xd/croot=1.24. When xd/croot>1.24, coefficient −CPV does not change with xd/croot anymore. xd/croot=1.24 is a proper place for the placement of TP, where the reversible energy transferring progress is completed and the calculated −CPV and CĖ do not contain the interference of the potential flow region. The remaining −CPV is about 13% of the total dissipation coefficient CΦ∗, which may indicate there are shock wave regions in the flow field. This term will be discussed next.Energy transfer occurs among decomposed components ofCĖ. The variation in the decomposed constituents of CĖ with xd/croot is shown in Figure 11. At the trailing edge of the aircraft (xd/croot=0), the absolute values of CĖa, CĖv, and CĖp are all over 50% of CD, while the overall CĖ sum is only 10% of CD. The signs of CĖa and CĖv are contrary to that of CĖp. Thus, the large absolute values of CĖa, CĖv, and CĖp are canceled in CĖ. The large absolute values of CĖa, CĖv, and CĖp near the aircraft trailing edge are caused by the potential flow region interference, which quickly decays within a fraction of croot.Figure 11
Sensitivity ofCĖa, CĖv, CĖp, and CĖ to the streamwise position of TP. The values of calculated dissipation components are normalized by near-field drag coefficient CD.Figures10 and 11 show the progress in energy transfer between the decomposed components of CΦ∗ and CĖ. Based on the variation charterer of −CPV and decomposed constitutions of CĖ, xd/croot=1.24 is a proper place for the placement of TP for the decomposition of CΦ∗, where the reversible energy transferring progress is completed and the calculated −CPV and CĖ components do not contain the interference of the potential flow region. At xd/croot=1.24, CĖa, and CĖp are nearly decaying to zero, and CĖ only accounts for the energy deposited in the trailing edge vortex, which is defined as the vortex dissipation component. All the viscous dissipation that occurred before xd/croot=1.24 was defined as the profile dissipation component. The no longer changing −CPV is treated as the volume pressure work-related dissipation components, which will be nonzero when the irreversible compression–expansion progress occurs [23].


## 4.2. Evaluation of BLI Benefits and Its Limitation of the BWB Aircraft with a BLI Propulsor

The power balance analysis framework presented in Section2.4. is used to quantify the BLI benefit of the BWB aircraft [27] and the D8 aircraft [21]. The results of the dissipation decomposition of the unpowered configuration are used in the analysis of the BLI benefits. The values of parameters fBLI, α, β, and ε, defined by Equations (24), (25), (27), and (29), are shown in Table 5, together with those of PSC, PSCmax, PSCjet, and PSCaf.Table 5
Comparison of PSC for the BWB aircraft and D8 [21].
fBLIαβεPSCPSCmaxPSCjetPSCafBWB22.5%9.0%60%0.35%18.8%3%2%D817.0%8.5%65%0.66%26.8%5%1%As shown in Table5, The power-saving coefficient PSC of the BWB aircraft [27] and D8 [21] is 5% and 6%, respectively. The evaluated limitations of the BLI benefit PSCmax of the two BLI aircrafts are 18.8% and 26.8%, respectively; this range is similar to the experiment results of Lv et al. [19, 20]. Despite having a higher value of fBLI, the value of PSC in the BWB aircraft [27] is still lower than that in D8 [21], which is mainly due to a lower value of ε.According to Equations (35) and (36), the term for the airframe dissipation reduction PSCaf and the term for the jet dissipation reduction PSCjet are largely determined by the propulsor design parameter ε. The value of PSCaf decreases with ε, while the value of PSCjet increases with ε. Thus, given the similar values of α and β, PSCjet is more dominant and PSCaf is weaker when the value of ε is higher. As indicated in Table 5, the value of ε for D8 is one-time greater than that for the BWB aircraft. As a result, the values of PSCjet and PSCaf are nearly the same on the BWB aircraft while in the PSC constitution of the D8 aircraft, PSCjet is the dominant factor. Additionally, the PSCjet of D8 is more evident than that of the BWB aircraft; that is, it is nearly one-time higher.


## 4.3. Energy-Based Aerodynamic Analysis of Aeropropulsion Coupling Effects

TheC2 case near the cruise point (CL=0.36, CX=0) is selected to conduct the dissipation decomposition. The cases of C0 and C1 with nearly a similar lift coefficient are chosen to conduct the dissipation constitution comparison with C2. In the CFD postprogress, PK is obtained from the mechanical flow power integral in Equation (3) on the inlet and outlet planes of the propulsion system, FX is obtained by integrating pressure and viscous force axial projection onto the airframe’s surface and axial momentum flux across the inlet and outlet planes, and lift force L is calculated by integrating pressure and viscous force vertical axial projection over the airframe’s surface. Total aerodynamic dissipation Φ∗ can be indirectly obtained through the power balance equation, Equation (8), or be calculated directly from the viscous dissipation Φ, volumetric pressure work PV, and the irreversible mechanical energy outflow Ė integral in Equation (7). The difference between the calculated Φ∗ from Equations (7) and (8) is within 5% of the value calculated from Equation (8) for all the three configurations. The vortex dissipation Φvortex in Equations (14) and ((17)) are calculated by the integral of Ė on the TP, according to Equation (17). The lift CL and net streamwise force CX coefficients are referenced to 1/2ρ∞V∞2Sref. The energy-related coefficients CPK, CΦ∗, CĖ, CPV, and CΦ are referenced to 1/2ρV∞3Sref, where Sref is the reference area.The results of the dissipation decomposition using Equation (7) are presented in Table 6. The ratio of CΦ to CΦ∗ of the three configurations is about 60%~70%, which contains the contribution from the aircraft surface boundary layer, trailing edge wake, and jet. The ratio of −CPV to CΦ∗ is quite considerable (10%~15%), which indicates there are shock wave regions in the flow field. The ratio of CĖ to CΦ∗ is about 18%~26%, lower than that of CΦ to CΦ∗.Table 6
Dissipation components forC0, C1, and C2.
CΦCĖ−CPVCΦ∗C00.0107(60.5%)0.0046(26.0%)0.0024(13.5%)0.0177C10.0190(66.7%)0.0052(18.2%)0.0043(15.1%)0.0285C20.0138(67.0%)0.0049(23.8%)0.(9.2%)0.0206To explain the variation inCΦ, the flow feature extraction methods described in Section 3.3. are used to isolate the jet region, the shock wave region, and the boundary layer/wake region in the flow field. The extraction priority is first the jet region, then the shock wave region, and last the shear layer region. The extracted regions are shown in Figure 12. The calculated viscous dissipation coefficients in the extracted regions are shown in Figure 13.Figure 12
Jet region, shock region, and shear layer region.Figure 13
Dissipation component comparison ofC0, C1, and C2.As is shown in Figure13, it is the difference in the dissipation that occurred in the shear layer region that causes the difference in CΦ. Then, the dissipation coefficient of shear layer CΦp is divided into five parts using the position-related sensor to show the variation in the dissipation generated on different parts of the BWB aircraft surfaces and in the wakes. The positions of the five dissipation components are labeled in Figure 14. Φcb‐bl represents the dissipation that occurred on the surfaces of the center body. Φow represents the dissipation that occurred on the surfaces of the blend wing and outer wing. Φeng‐bl represents the dissipation that occurred on the surfaces of the propulsor nacelle. Φcb‐wake represents the dissipation that occurred in the center body trailing edge wake. Φeng‐wake represents the dissipation that occurred in the propulsor nacelle trailing edge wake.Figure 14
The position ofΦcb‐bl, Φow, Φeng‐bl, Φcb‐wake, and Φeng‐wake.
(a)(b)According to the results of the decomposition of the shear layer dissipation coefficientCΦp of C0, C1, and C2 in Table 7, the dissipation components with obvious variation are CΦcb‐wake, CΦeng‐bl, and CΦeng‐wake. By comparing the CΦcb‐wake term of C1 and C2, the wake dissipation reduction brought by BLI is 0.0005, while that calculated from Equations (36) and (32) is 0.0004. These results show that the airframe dissipation reduction is underestimated by Equation (36), mainly due to the fact that we ignore the change in the velocity profile of the trailing edge of the center body due to propulsor installation. Meanwhile, the variations in CΦcb‐bl and CΦow of the three configurations are less obvious. The value of CΦcb‐bl of C2 is slightly less than that of C0 and C1 due to the reduction in the surface area of the center body. The large increase in CΦeng‐bl and CΦeng‐wake may be due to the flow separation on the propulsor outer surface as is shown in Figure 15, which can be reduced through a parametric design of the configuration.Table 7
Shear layer dissipation coefficientCΦp decomposition.
CΦcb‐blCΦcb‐wakeCΦeng‐blCΦeng‐wakeCΦowCΦpC00.00520.0004000.00450.0101C10.00520.00100.00330.00160.00440.0155C20.00480.00050.00080.00020.00440.0107Figure 15
The contours of Mach number and streamline ofC1.To explain the variation inCΦcb‐wake and CΦeng‐wake, the trailing edge velocity profiles of the three configurations (C0, C1, and C2) have been extracted from the previously obtained CFD results, which are presented in Figure 16. The extraction positions of the trailing edge velocity profiles are marked in Figures 16(a)–16(c). Figures 16(d)–16(f) are the trailing edge velocity profiles of the center body and propulsor for C0, C1, and C2, respectively. The center body trailing edge velocity profile of C1 (Figure 16(e)) is steeper than that of C0 (Figure 16(d)), which represents a larger velocity gradient of the boundary layer flow at the trailing edge and thus a larger CΦcb‐wake. The center body trailing edge velocity profile of C2 (Figure 16(f)) is moderate compared with that of C1 (Figure 16(e)), which is the result of the wake filling effect of BLI. The moderate center body trailing edge velocity profile of C2 represents a smaller velocity gradient of the boundary layer flow and brings a reduction in CΦcb‐wake. The reduction in CΦeng‐wake of C2 compared with C1 can also be attributed to a moderated trailing edge velocity profile.Figure 16
Trailing edge velocity profiles forC0, C1, and C2. (a–c) Positions of the extracted trailing edge velocity profiles for C0, C1, and C2. (d–f) Trailing edge velocity profiles of center body and propulsor for C0, C1, and C2.
(a)(b)(c)(d)(e)(f)The volumetric pressure work coefficient−CPV of C0, C1, and C2 is 0.0024, 0.0041, and 0.0019, respectively. The value of −CPV of C2 is smaller than that of C1, which means a reduction in the waste of mechanical energy. As is explained by Sanders and Laskaridis [23], the value of −CPV represents the mechanical energy waste in the control volume that results in an increase in internal energy. This term cannot be neglected when there is a shock wave in the flow field. When the fluid passes through a shock wave region, the mechanical energy will decrease because of expansion damping [23].In order to show the influence of the shock wave on the variation of−CPV, the flow fields at Ma=0.85 and Ma=0.6 of C0 were generated by the CFD method. Then, the 2D slices of the two flow fields at the spanwise location Z=2m were extracted, and the contours of velocity divergence ∇⋅V and relative pressure contours p−p∞ are shown in Figure 17. In this figure, the contour flood is for ∇⋅V, and the contour line is for p−p∞. Figure 17(a) is extracted at Ma=0.85, and Figure 17(b) is extracted at Ma=0.6. As is shown in Figure 17(a), there is a strong compression region in the high Mach number (Ma=0.85) flow, identified by the large negative value of ∇⋅V. This strong compression region is actually the shock wave region, and there may be obvious mechanical energy loss in this region. Thus, the contribution of −CPV to CΦ∗ cannot be neglected. However, in Figure 17(b), there is no obvious compression region in the flow field with Ma=0.6. Thus, the contribution of −CPV to CΦ∗ is only 4%, which can be neglected. The above comparison shows that −CPV is highly affected by the existence of the shock wave region and, thus, is highly related to the mechanical energy loss due to the shock wave. The obvious variation in −CPV among C0, C1, and C2 may due to the change of strength of the shock wave.Figure 17
The contour of velocity divergence∇⋅V and relative pressure p−p∞. The contour flood is for ∇⋅V, and the contour line is for p−p∞. The labeled figures are the value of p−p∞. (a) The contours at Ma=0.85. (b) The contours at Ma=0.6.
(a)(b)The coefficients of the dissipations that occurred in the trailing edge vortexCĖ of C0, C1, and C2 are 0.0046, 0.0052, and 0.0049, respectively. The variation in CĖ between C1 and C2 shows that the installation height of the propulsor system can change the lift coefficient distribution along the aircraft spanwise, which is the determining term of lift-induced vortex strength. To explain the variation in CĖ of C0, C1, and C2, the contours of Ev_uni on TP located at X=100m (2.3croot) of C0, C1, and C2 are generated from the CFD results, as shown in Figure 18. Ev_uni is the nondimensionalized transverse kinetic energy deposited in the trailing edge vortex, and it is defined as
(42)Ev_uni=1/2v2+w21/2V∞2.Figure 18
Ev_uni contours of C0, C1, and C2. (a–c) TP plane positions. (d–f) Ev_uni contours on TP plane.
(a)(b)(c)(d)(e)(f)The dark red region in Figure18(e) represents a notable contribution to Φvortex of C1, which is weaker in the Ev_uni contour of C0 and C2 (Figures 18(d)–18(f)). The difference in Ev_uni distribution of C0, C1, and C2 shows the source of variation in Φvortex in all three configurations.Summarizing the abovementioned analysis on the results of the total dissipation coefficientCΦ∗ breakdown of C0, C1, and C2 for the BWB aircraft shows that the airframe dissipation components that are strongly affected by BLI effect are −CPV, CΦwake, and CΦvortex. Variations in CΦeng‐bl and CΦeng‐wake are not strictly caused by BLI effects, but they are a result of inappropriate design of the propulsor nacelle that can be improved later.


## 5. Conclusions

An energy-based aerodynamic analysis method is used to assess the BLI benefits applied on a BWB aircraft. The power-saving coefficient PSC is chosen as the metric for the calculation of BLI benefits. The determining expression of PSC is further related to the BLI fractionfBLI, two airframe dissipation ratios α and β, and the propulsor specify thrust ε. By letting fBLI=1, the maximum of power-saving PSCmax can be obtained at fixed α, β, and ε. PSC and PSCmax will increase with higher ε and fBLI. In addition, the contribution of jet loss reduction in PSC is more dominant with a higher ε. The evaluated PSC of the BWB aircraft with a BLI propulsor is 5% using the previously developed determining expression of PSC. The maximum PSC is 18.8%.The CFD calculations of one power-off and two power-on configurations are conducted to obtain the aerodynamic dissipation components of the aircraft in detail. The three configurations are namedC0, C1, and C2. C0 is the unpowered configuration. C1 is the power-on configuration with a non-BLI propulsor. C2 is the power-on configuration with a BLI propulsor. The total aerodynamic dissipation of the aeropropulsive coupling flow around the aircraft is broken down into three viscous dissipation components and one volumetric pressure work. All of the three viscous dissipation components are related to a specific flow feature. The main findings of the analysis of dissipation components are listed as follows:
(1)
The BLI benefit evaluation based on the dissipation components of the unpowered configurationC0 slightly underestimates the power saving of BLI. The reason is that the variation in the boundary layer velocity profile at the trailing edge of the center body part is neglected(2)
Φvortex has been changed by the BLI effects given that the center body part also provides lift for the aircraft(3)
The volumetric pressure work term−PV is strongly affected by the BLI effects, and it is neglected in the previous BLI analysis and should be analyzed comprehensively in the future



---

*Source: 6452099-2022-06-16.xml*