Computational fluid dynamics challenges for hybrid air vehicle applications

This paper begins by comparing turbulence models for the prediction of hybrid air vehicle (HAV) flows. A 6 : 1 prolate spheroid is employed for validation of the computational fluid dynamics (CFD) method. An analysis of turbulent quantities is presented and the Shear Stress Transport (SST) k-ω model is compared against a k-ω Explicit Algebraic Stress model (EASM) within the unsteady Reynolds-Averaged Navier-Stokes (RANS) framework. Further comparisons involve Scale Adaptative Simulation models and a local transition transport model. The results show that the flow around the vehicle at low pitch angles is sensitive to transition effects. At high pitch angles, the vortices generated on the suction side provide substantial lift augmentation and are better resolved by EASMs. The validated CFD method is employed for the flow around a shape similar to the Airlander aircraft of Hybrid Air Vehicles Ltd. The sensitivity of the transition location to the Reynolds number is demonstrated and the role of each vehicle£s component is analyzed. It was found that the ¦ns contributed the most to increase the lift and drag.


INTRODUCTION
Predicting the vortical §ows present around lighter than air (LTA) and HAVs [1] is challenging for CFD.In the literature, prolate spheroids have been commonly used to approximate LTAs [2,3], since conventional airships have ellipsoidal axisymmetric shapes [4].The present work suggests that this is true; there are, however, substantial di¨erences.When spheroids and airships are pitched at high angles, vortices are shed due to their shape curvature and pressure gradients.In the case of spheroids, these vortices are very close to each other leading to vortical interactions further downstream.Hybrid air vehicles, on the other hand, tend to be much wider than spheroids and, especially in modern HAV designs, the hull is composed by more than one lobes to maximize the production of lift as shown in Fig. 1.This makes the distance between vortices much longer than for prolate spheroids and leads to a much weaker or even nonexistent interaction between them.In addition, airships have ¦ns and/or strakes that ¦x the point of vortex formation and enhance their strength.
Computational §uid dynamics is typically based on one-or two-equation turbulence models and a Boussinesq£s approximation [5].These models assume an alignment between the Newtonian and Reynolds stress tensors.For many aerospace applications, this assumption is valid.However, for more complex problems, such as spheroids or airships, where the §ows are highly threedimensional (3D) and have strong secondary §ows, this assumption is very restrictive.Boussinesq models also tend to over-produce eddy viscosity near vortex cores, unless some clipping is introduced in their production terms.For these cases, more advanced turbulent models are required.One option is the use of Reynolds-stress Models (RSMs), where a transport equation is solved for each component of the Reynolds stress tensor [5].The implementation of such models is complex and their use is costly for routine applications.Few works in the literature are found where RSMs are employed for spheroids.Alpman and Long [6] used the PUMA solver for the analysis of the §ow around a 6 : 1 spheroid and a sphere.The results showed good agreement with experiments with some discrepancies in the skin friction predictions.Although the computations were costly in terms of CPU resources, the analysis of the normalized turbulent stresses in the sphere demonstrated its highly anisotropic nature and, therefore, justi¦ed the use of RSM.More recently, Yossef et al. [7] studied the spheroid at 20 • of pitch and 4.2 • 10 6 Reynolds.Boundary layer velocity components and Reynolds stresses obtained with RSM showed good agreement with the experiments [8].
Explicit algebraic RSMs (EARSMs) [911] are the nonlinear models that lay in between Boussinesq models and RSM.Explicit algebraic RSMs are based on two-equation turbulence models, but replace Boussinesq£s relationship between the Reynolds stress tensor and the main strain rate with a higher-order expansion.The isotropy restriction of the Boussinesq models is, therefore, lifted.In the literature, Grudestam et al. [10] employed an EARSM for the analysis of rotating two-dimensional (2D) shear and channel §ows and 3D pipe §ows.Boundary layers under adverse pressure gradient, wakes, and mixing layers were studied by Hellsten [11].
There are experimental data available in the open literature for 6 : 1 prolate spheroids [8,12,13] and for 3 : 1 spheroids [14].Scott [15,16] studied the 6 : 1 spheroid [12] at 20 • of pitch using the OVERFLOW solver.The ¦rst study [15] showed the importance of grid re¦nement in capturing the vortical structures.The kω SST model showed better pressure agreement with experiments than the kω and SpalartAllmaras (SA) models.However, the skin friction was underpredicted by all the models.Detached Eddy Simulation (DES) was employed in a later study [16]; however, no improvements in the pressure and skin friction predictions were reported.This was attributed to the lack of separation of the main vortices at 20 • of pitch, with the DES model working as Unsteady RANS (URANS) only.
Similar DES results were obtained by Xiao et al. [17].Although primary and secondary separations were in agreement with the experiments [12], there was no improvement in the skin friction predictions.Similar studies were performed by other authors [1820].
Sorensen [21] also employed the 6 : 1 spheroid [12] for the validation of the γRe θ transition model, for Reynolds between 3.2 • 10 6 and 9.6 • 10 6 , and at 0 • and 30 • of pitch.Fully turbulent and transitional cases showed good agreement with experiments for the pressure.At zero pitch, the transition model agreed better with the experimental skin friction than the fully turbulent cases, but the γRe θ model was not able to predict the correct location of transition at 30 • .This was attributed to the lack of cross- §ow transition prediction capability of the employed model.
Vizhinho et al. [22] developed a transition model based on SA and pretransitional quantities (V-SA).The experiments were compared against the baseline SA, γRe θ model, and the developed V-SA transition model.The latter predicted the closest transition onset to the experiments.However, the predicted transition length was shorter than in the experiments and was attributed to an excess of turbulent kinetic energy di¨usion inside the boundary layer.
The explicit algebraic model developed by Rumsay and Gatski [23] was employed by Morrison et al. [24] for the analysis of a 6 : 1 spheroid at 30 • of pitch and at Reynolds of 43 • 10 6 .Better agreement in the pressure predictions was obtained with the EASM than with the kω model [25].However, lack of grid convergence was reported and skin friction predictions were not investigated.This paper explores di¨erent turbulence models aiming to ¦nd the best framework to model airship-like §ows.A study of the aerodynamics for a variant of a HAV is then provided.

HMB2 Flow Solver
The Helicopter Multi-Block (HMB2) CFD solver [26] is used for the present work and has so far been validated for a number of applications, including helicopters, wind turbines, turboprop, and UCAV (unmanned combat aerial vehicle) aircraft.HMB2 solves the NavierStokes equations in integral form using the Arbitrary Lagrangian Eulerian formulation for time-dependent domains with moving boundaries: where V (t) is the time-dependent control volume; ∂V (t) is its boundary; w is the vector of conserved variables [ρ, ρu, ρv, ρw, ρE] T ; and F i and F v are the inviscid and viscous §uxes, including the e¨ects of the mesh movement.The Navier Stokes equations are discretised using a cell-centered ¦nite-volume approach, leading to the following equation: where i, j, and k are the cell indices; w are the variables; R are the residuals; and V is the volume.The upwind scheme [27] is used for the discretization of the convective terms and MUSCL (monotone upstream scheme for conservation laws) [28] variable extrapolation is used to provide formally 3rd-order accuracy.
To account for low-speed §ows, the low-Mach Roe scheme of Rieper [29] is employed [30].The linearized system is solved using the generalised conjugate gradient with a block incomplete lower-upper preconditioner [31].

Turbulence Modeling
Fully turbulent §ows are usually assumed in aerospace applications; however, some problems involve both laminar and turbulent §ows that need to be correctly captured.For this type of problems where the e¨ect of transition is important, transition models such as γRe θt should be employed.On the other hand, EASMs (EARSMs) provide good prediction of the vortex cores.In addition, unsteady problems with large separation can be predicted with Scale-Adaptive Simulation (SAS).The HMB2 solver has a library of turbulence closures which includes several one-and two-equation turbulence models and versions of the kω model, including the SAS model [32], the kω EASM (EARSM) [11], and the γRe θt transition model [33].

Explicit algebraic stress kω model
The following are the equations for the kω baseline (BSL) model of Menter [34].
In the EARSM model from Hellsten [11], the production term (P ) and the eddy viscosity (µ t ) have additional nonlinear contributions: The terms a¨ected in the production (P ) are the viscous stresses (τ ij ), as presented in Here, where One contribution comes from the addition of the anisotropic terms (a The other contribution comes from the scaling of the eddy viscosity (µ t ) by a nonlinear coe©cient C µ , given in In this model, β * = 0.09, like in the BSL model [34].The closure coe©cients γ, β, σ k , σ ω , and σ d are obtained using the following equation with the blending function f b : Here, where Note that the terms • 1 and • 2 are exactly the same as the ones in the BSL model.The coe©cients for the blending function are detailed in [11].

Scale-adaptive simulation
Scale-adaptive simulation [32] is a modi¦cation of an SST-RANS model [34] based on the use of the second mechanical scale (in the form of second derivatives of velocity) in the source terms of the underlying turbulence model.An additional term Q SAS is then added in the transport equation for the speci¦c dissipation rate ω, where l is the length scale; l vK is the von Karman length scale in classic boundary layer de¦nition; and the constants ζ 2 = 3.51, σ = 2/3, and C = 2. Scale-adaptive simulation is a hybrid RANS / Large-Eddy Simulation (LES) model which can produce spectral content for unstable §ows, adjusting the turbulence length scale to the local §ow inhomogeneities and balancing the contributions of modeled and resolved parts of the turbulent stresses.For steady §ows, it acts as a RANS model and for §ows with transient instabilities like those with massive separation, the model reduces its eddy viscosity according to the locally resolved vortex size represented by the von Karman length scale.The SAS model can resolve the turbulent spectrum down to the grid limit and avoids RANS-typical single-mode vortex structure.

Local correlation-based transition model
The γRe θ , or Local Correlation-based Transition Model (LCTM) [35], uses two additional transport equations: one for intermittency (γ) and one for the transition momentum thickness Reynolds number (Re θt ), which is formulated in terms of the scalar quantity, Re θt .The γ-equation is used to trigger the transition process and the Re θt -equation is employed to avoid nonlocal variables and pass information from the freestream to the boundary layer.The transport equations are: where ρ is the density; u j is the velocity vector; µ is the molecular viscosity; and µ t is the eddy viscosity.The production terms, P γ1 , P γ2 , and P θt are de¦ned as follows: where Ÿ is the vorticity magnitude; S is the strain rate magnitude; U is the local velocity magnitude; the parameters F length and F onset are used to control the length and onset location of transition, respectively; and F turb and F θt are the parameters for controlling the destruction/relaminarization of the boundary layer and the boundary layer detector, respectively.Term P γ2 acts as a sink and ensures that the intermittency remains close to zero in the laminar boundary layer.The production term P θt in Eq. ( 3) is designed to force the transported scalar Re θt in the freestream to match the local value of Re θt into the boundary layer.Parameter Re θt is calculated locally from the empirical correlation proposed by Menter et al. [35].This parameter is the critical Reynolds number where the intermittency ¦rst starts to increase in the boundary layer.This occurs upstream of the transition Reynolds number.Likewise, Re θt is the location where the velocity pro¦le begins to deviate from the purely laminar pro¦le.The boundary condition for the intermittency factor (γ) at the freestream and wall boundaries is zero §ux.For " Re θt , the boundary condition at the freestream and wall boundaries is also zero §ux.Finally, the model constants are c a1 = 2; c e1 = 1; c a2 = 0.06; c e2 = 50; c a = 0.5; σ γ = 1; σ θt = 2; and c θt = 0.03.

MESH GENERATION
Multiblock structured topologies were employed to allow for good representation of the spheroid and airship surfaces.The blocks are also used for easy sharing of the computational load between parallel processors.
As shown in Fig. 2, an O-topology was employed for the 6 : 1 prolate spheroid grid, with 576 blocks for good load balancing between processors.The in §ow, out §ow, and far-¦eld boundaries extended 10L from the spheroid surface where L is the length of the spheroid.Since the spheroid is symmetric with respect to its long axis, only half of it was meshed.A ¦ne grid was employed for good boundary layer resolution and for capturing the vortices found above the spheroid.Hence, the o¨-body grid was kept uniform in the region where vortices were expected, with maximum cell sizes of 0.03%L (Fig. 2b).Good grid resolution close to the wall is also a requirement for the γRe θ turbulence model.In the ¦rst block normal to the surface, 74 cells were used, with a ¦rst cell spacing of 10 −5 L and an expansion ratio of 1.05, to ensure y + < 1. Around the spheroid azimuth, 480 cells were employed and 592 cells were used along its long axis, with a maximum cell  size in that direction of 0.003L.This led to a 40-million cells grid.Assuming symmetry, only half of the HAV£s hull was meshed.At the in §ow, out §ow, and far-¦eld boundaries, freestream conditions were imposed.Five sliding planes [36] were employed to allow for localized grid re¦nement near the body and to have a cartesian grid in the rest of the domain.These planes are the ones close to the body shown in Fig. 2c.An O-topology was employed around the HAV for optimal orthogonality to the surface (Fig. 2e).The ¦rst cell size of 3 µm was employed to ensure y + < 1 and 50 cells were used to capture the boundary layer.Around the hull, 516 cells were employed and 362 cells along its length.The grid had a total of 31 • 10 6 cells.
Four con¦gurations were considered to study the role of the components of an approximated hybrid air vehicle (Fig. 3).Con¦guration AL 1 consists of the bare hull; the ¦ns are then added (AL 2) and also the leading edge extension (LERX) in con¦guration AL 3. Finally, strakes are considered in con¦guration AL 4. To assess grid convergence, two grid levels were employed for the baseline case (grids AL 4 and AL 5).

RESULTS AND DISCUSSION
In this section, the computational results are presented.For the 6:1 spheroid computations, a Reynolds number of 4.2 • 10 6 (based on the spheroid£s length L) and Mach number of 0.15 were employed.In addition, the spheroid was pitched at 20 • .For the analysis of the HAV hull, a Reynolds number of 3 • 10 6 (based on the vehicle£s length) and a Mach of 0.12 were selected.The HAV was pitched at 20 • and 30 • .Further computations are then presented for the approximate HAV (AL), including the role of the strakes in the aerodynamic coe©cients and e¨ect of the §ow conditions and Reynolds number on the transition onset.For these cases, a wind speed of 40 m/s was employed and sea level conditions were assumed.Unless otherwise speci¦ed, the Reynolds number, based on the vehicle£s length (L), was 3 • 10 6 .The RANS and URANS computations were performed, with the kω SST turbulence model by Menter [34].Table 1 provides a summary of the test cases, including the employed grids, turbulence model, and the number of CPUs used in the computations.The ¦rst computations were performed for the study of di¨erent turbulence modeling frameworks.The last 5 sets of computations were performed for the analysis of the components of the HAV.
The spheroid §ow at high pitch is characterized by the development of vortices due to curvature e¨ects.Figure 4 presents a schematic of the §ow direction and  1) is created at an azimuth of ϕ = 160 • .Secondary (2) and tertiary (3) vortices also appear in this type of §ow, located at ϕ = 140 • and 130 • , respectively.Indication of their rotation direction (as seen from the front looking downstream) is included.

Validation of the Method
The HMB2 method is ¦rst validated with the pressure measurements performed by Wetzel [12] and the friction measurements by Chesnakas et al. [13].In the experiments [12], the spheroid was tripped at 20%L.Pressure taps were installed and the model was rolled to sweep the transducers from the windward to the leerward sides and to map out the pressure distribution over the model surface.Data were taken only for one half of the model, assuming symmetry of the §ow.The skin friction data were not measured directly; the velocity pro¦les were measured and ¦tted to a Spalding-type wall law [13].
Figure 5 shows pressure and friction coe©cients around the spheroid azimuth, through a slice at x = 77% L. The CFD results (kω SST [34], kω BSL + EARSM [11], and kω-γ-Re θ transition model [35]) and experiments [12,13] are compared.The DES results presented by Constantinescu et al. [19] are also included.The SST computation was performed using 40-and 12-million cells grids.The solutions on both grids are practically identical and, therefore, grid independence was assessed in terms of pressure and friction coe©cients.For the rest of the paper, however, the ¦ne grid was employed as it enabled better resolution of the vortices o¨-body.
The presence of the main vortex (V1) at ϕ = 160 • very close to the spheroid surface leads to a drop in pressure as can be observed in Fig. 5a, which is in very good agreement with experimental data [12].A small drop in pressure is also observed at ϕ = 140 • , corresponding to the secondary vortex (V2), and a small peak at 130 • of azimuth, where the tertiary vortex (V3) is located.The di¨erences between the SST and the EARSM models are very small in the lower surface of the spheroid (from 0 • to 90 • ).Once the vortices are developed, the EARSM improves the agreement with experiments.The transitional model, however, presents poor agreement with experiments.This is expected, since, during the experiments, the spheroid was tripped at 20% L and at this location, the boundary layer was fully turbulent.The kω-γRe θ model predicts free transition and at this axial location, parts of the boundary layer are laminar.The behavior of this model is further discussed in subsection 4.5.
For the friction coe©cient (C f ) presented in Fig. 5b, a much better agreement with the experiments [13] was obtained with the EARSM, mainly in the region where the secondary (V2) and tertiary (V3) vortices are created at ϕ = 140 • and 130 • , respectively.Note that this prediction is much closer to the experiments than the one reported by Constantinescu et al. [19].At azimuth of 160 • , where the main vortex (V1) is located, all models underpredict the peak in C f .This peak underprediction was also reported in [15,16,19], even though they employed di¨erent turbulence models (RANS and DES with the SA model).These discrepancies in friction could be due to di©culties of the turbulence models to correctly predict this vortical structure, which seem to be weaker in the CFD, but also di©culties in the experiment to measure these quantities in the wind tunnel with high level of accuracy.Computations on a full spheroid without symmetry condition did not improve this aspect of the comparison.Much the same way, di¨erent laws for evaluating the C f including ¦ts to velocity with polynomials showed little sensitivity due to the employed ¦ne mesh.

Analysis of the Stress Tensor
To better understand the mechanism of the EARSM, this subsection provides an overview of the components of the turbulent stress tensor and a comparison with the kω SST model.As presented in subsection 4.1, the EARSM improved the agreement with the experiments with respect to kω SST, due to a reduction of the eddy viscosity in the vortex core.This can be seen in Fig. 6 that compares the turbulent eddy viscosity (Re T = µ T /µ) through a slice at x = 0.77 % L. As can be seen, at the edge of the primary vortex, the eddy viscosity levels are approximately the same for both turbulence models.On the other hand, a strong reduction in the eddy viscosity at the core is observed in the EARSM (Fig. 6c).This reduction of the turbulence levels is due to the nonlinear factor (C µ /β * ) that scales the turbulent eddy viscosity, depending on the local shear and vorticity, as was shown in Eq. ( 2).For the same axial slice through the main vortex, Fig. 6d shows the distribution of the C µ /β * scaling factor where the eddy viscosity at the core is completely eliminated.
In the EARSM, the stress tensor is de¦ned as where τ * ij is the Boussinesq stress tensor but employing the nonlinear turbulent viscosity shown in Eq. ( 1) and a (ex) ij is the anisotropy tensor.This is a symmetric matrix, whose components are shown in Fig. 7, with contours on a slice through the main vortex (at x = 77% L).Note that a (ex) ij is scaled with 1/(β * ω). Figure 8 shows the ¦rst three invariants of the two contributions of the total stress tensor of Eq. ( 4): the stress contribution (τ * ij ) and the anisotropy contribution (a (ex) ij ρk).As can be seen in Figs.8a and 8d, the ¦rst invariant of the anisotropy tensor (I a = tr(a)) is practically zero and, therefore, the main contribution to the stress tensor comes from τ * ij .The second, II a = (1/2) tr(a) 2 − tr(a • a) , and the third, III a = det(a), invariants show that the anisotropy is normal to the viscous stress as for the same cell, each contribution is of the opposite sign.
A comparison of the cross terms of the stress tensor between the EARSM and kω SST models is presented in Fig. 9.The values were extracted along two circular lines: one passing through the main vortex (see Fig. 9, left column) and another line through the core of the secondary vortex (see Fig. 9, right column).As can be observed, the peaks are located at the same azimuthal position for both models (ϕ = 140 • for the secondary vortex and ϕ = 160 • for the main vortex).In general, larger amplitudes are observed in the EARSM that seems to capture more spatial variations than the kω SST.This can be seen in Fig. 9, right column.

Comparison between URANS-EARSM and SAS-SST models
The URANS results employing the SST model [34] and EARSM [11] showed that the later was able to capture better the vortices generated around spheroids, leading to a better agreement with the experimental data.This subsection aims to make a comparison between the URANS-EARSM and SAS-SST [32] computations.
Figure 10 shows the slices of pressure and friction coe©cients at x = 0.77% L. For the SAS-SST model, two time-step sizes were employed (-t * = 0.02 and 0.005) and identical solutions were obtained, showing independence on the size of the time step.The SAS-SST model did not improve the agreement with the experiments unlike the EARSM model.This is due to the fact that the SAS model reduced the overall eddy viscosity on the vortex but not locally at the core.This can be observed if the turbulence levels of the two models are compared (see Figs. 6b and 11a).A line crossing the main vortex (Fig. 11b) shows that the eddy viscosity at the edge of the vortex is higher for the EARSM model which leads to stronger vortices.These results show that to predict this §ow, accounting for anisotropy is essential.

Analysis of the Flow
Figure 12 shows contours of Q-criterion [37] at di¨erent spanwise sections along the spheroid and a detail at x = 77% L. Friction over the spheroid is also included.As can be observed, a main vortex is generated located at 160 • of azimuth.The secondary vortex is also captured, located at an azimuth of 140 • , which is in good agreement with experiments [13].The third vortex is also present at 130 • , but is much weaker than the other two.As the EARSM model produces less eddy viscosity at the vortex core, the vortices are tighter and less di¨used at their edge than with kω SST or the SAS-SST models.This will lead to a stronger interaction of vortices at the rear of the spheroid.The SAS-SST model presents lower levels of Q-criterion as a result of the overall reduction of eddy viscosity throughout the §ow-¦eld.
In Fig. 13, the contours of helicity (h) on a slice at x = 77% L are presented.Helicity is obtained by the dot product of the velocity and the vorticity vectors, h = U U U • (∇ × U U U ). Good agreement with experiments [13] is observed.As can be seen, the secondary vortex is further detached from the wall when the EARSM is employed.Likewise, lower levels of helicity in the secondary vortex are predicted by the SAS-SST model.
Friction lines projected onto a plane are shown in Fig. 14.A close view zoom of a region on the upper surface (from 90 • to 180 • of azimuth) is also included that extends from mid-length to the rear of the spheroid.As can be observed, the primary separation line (A) is practically identical in both URANS solutions (ϕ = 105 • at 77% L).Conversely, the SAS-SST model predicts the primary separation at higher azimuth angle (ϕ = 105 • at 77% L).There are di¨erences between all three models in the secondary separation line (B).
The EARSM predicts a line that begins at 55% L and stays at 150 • .The SST model predicts a later separation (from 65% L) and moves the separation line inboards from 140 • to 150 • further downstream.For the SAS model, it also starts from 65% L but in this case, stays at approximately 150 • .This earlier secondary separation predicted by the EARSM also reported by Morrison et al. [24].Very small di¨erences are observed for the reattachment line (C) between the URANS solutions, while the SAS model predicts a line at later azimuth.
In Fig. 15, the primary and secondary separation lines location is validated against experiments [13].These lines were measured using constant temperature and constant current anemometry and using oil §ow visualization.It should In the present CFD results, a fully turbulent §ow was assumed.Di¨erences between experiments and CFD from 0 to 30% L are expected and results should be only compared from 30%L.Similar trends between CFD and experiments are observed and the best agreement for both primary and secondary lines is obtained between the EARSM and the oil §ow measurements.

E¨ect of Transition
In this subsection, transition e¨ects for the spheroid pitched at 20 • are studied with the kωγRe θ transition model of Menter [35].  1, respectively).For the later, the solutions at three instances within one §ow travel time are provided.Instance t 1 corresponds to a phase of ψ = 0 • of the oscillation of the transition boundary, instance t 2 corresponds to a phase of 120 • , and instance t 3 corresponds to ψ = 240 • .The case where the SAS-SST model was used in Fig. 16a clearly shows a fully turbulent §ow.As can be seen in Fig. 16b, most of the lower surface of the spheroid remains laminar if the γRe θ is employed.In addition, due to the large pitch angle and the surface curvature, the turbulent region moves from an azimuth angle of 180 • at 5% L to 95 • at 50% L and 60 • at 95% L. However, this transition boundary does not stay ¦xed with time as shown in Fig. 16b.It should be noted that unlike the previous cases, the separation and reattachment lines are also not ¦xed with time.
Figure 17 presents surface friction lines at three instances within one §ow travel time.As can be observed, the primary separation line (A) stays practically the same from nose to 80% L. However, at instance t 3 , the streamlines change curvature and the separation line from 80% L to the rear change azimuth angle from 80 • to 70 • .The secondary separation line (B) stays around 150 • , but oscillates with ±5 • of amplitude.Compared to the fully turbulent cases (see Fig. 14), this secondary separation line is located much further downstream.Finally, the reattachment line (C) is similar at all instances.
This unsteadiness a¨ects the location of the vortices.For the same instance within one travel time (t 3 ), Fig. 18a shows the variations in C f on a slice through x = 77% L. Azimuth angles of 0 • to 90 • correspond to the lower surface and the boundary layer was predicted as laminar.In this region, the skin friction remained constant.The fully turbulent region (from 150 • to 180 • ) also stayed constant through time.Oscillations are observed in the region in-between (from 90 • to 150 • ) due to the oscillations of the transition location with time.This is the result of the expansion and contraction of the main vortex.As can be seen in Fig. 18b, the main and secondary vortices moved closer to the spheroid surface at instance t 3 .These results show the importance of ¦xing the transition point to fully characterize the §ow around this type of shapes.

E¨ects on the Loads
Table 2 provides a summary of the drag and lift for the studied cases.Note that the coe©cients are nondimensionalized taking the URANS kω SST as reference.Since an oscillatory pattern was observed in the transitional case, the loads presented for this case were averaged.Taking the URANS kω SST as baseline, the total drag predicted by the EARSM increased by 3% which was due to an increase in the pressure drag, while the friction drag stayed practically the same.Using the SAS-SST formulation, the total drag was reduced by 7.5%.However, the predicted friction drag in this case was higher and the pressure drag lower than for the URANS solution.The transition model predicted a much lower overall drag (−34%), where both the C dp and C d f were reduced.This These results show that although there is substantial di¨erence in the vortical structures predictions and how these a¨ect the pressure and friction on the spheroid£s surface, the overall loads are not a¨ected as much for the fully turbulent cases.The e¨ect of transition conversely is very important and highly a¨ects both lift and drag predictions.

Application Case ¡ Bare Hull of a Hybrid Air Vehicle
The analysis on spheroids showed the di¨erences between turbulence models.As an application example, the bare hull of a HAV approximating the Airlander of Hybrid Air Vehicles Ltd. at model scale (see Fig. 1b) is studied here.

Analysis of turbulent quantities
Figure 19 shows the contours of turbulent eddy viscosity on a slice at 80% L. As can be observed, at 20 • of pitch, the EARSM slightly reduces the turbulence levels.This is better observed in Fig. 20a that shows a comparison of Re T on lines passing trough the main vortex.Note that the shaded vortices included in this ¦gure are the contour lines of Q-criterion.These small di¨erences between models are due to the fact that at this pitch angle, the vortices are not fully rolled up.If the pitch angle is increased to 30 • , the vortices are fully developed and the EARSM behaves similarly to the spheroid cases.The eddy viscosity at the vortex core is in this case highly reduced, as can be seen in the contours in Fig. 19b (right column) and in the lines through the main vortex in Fig. 20b.

Analysis of the §ow-¦eld
Slices with contours of Q-criterion along the HAV are provided in Fig. 21, with a detail of the §ow at 80% L. At this axial station, only the main and sec-   Friction lines projected onto a 2D plane are presented in Fig. 22.Only one separation line (A) on the upper surface is observed at 20 • of pitch starting Figure 22 Friction lines projected onto a 2D plane at Re = 3 • 10 6 , M = 0.12, and α = 20 • (a) and 30 • (b): 1 ¡ kω BSL + EARSM; and 2 ¡ kω SST at 70% L, with the same pattern for both turbulence models.This di¨ers from the spheroid §ow that at this pitch angle had primary and secondary separation lines and a reattachment region, as seen in Fig. 14.At 30 • of pitch shown in Fig. 22b, a similar pattern to the spheroid£s is observed.There is a 3 • of azimuth di¨erence between turbulence models in the location of the primary separation line (A).The secondary separation (B) and reattachment (C) lines are predicted by both models at the same axial and azimuthal positions.
Figure 23 compares the primary and secondary separation lines between the spheroid and the HAV at α = 20 • and 30 • .There is very little resemblance between spheroid and HAV at 20 • of pitch.In the latter, only primary separation is present and starts much further downstream than for the spheroid.On the other hand, a very similar pattern is observed between spheroid and HAV at 30 • of pitch.It is interesting that both primary and secondary separation lines are

Analysis of the aerodynamic loads
Figure 24 shows azimuthal surface pressure and friction distributions on a slice at 80% L. Note that the analysis is focused on the outer lobe of the HAV where the vortices are developed.The layout of the vortices with isolines of Q-criterion is included in the ¦gures.Azimuth ψ = 0 • corresponds to the lower surface, ψ = 90 • is on the side, and ψ = 180 • is the uppermost point in the upper surface.At 20 • of pitch, as shown in Fig. 24a, the EARSM model presents a dip in pressure and a high peak in friction at ϕ = 180 • which is an indicative of the generation of a vortex.For this case, the results were compared with a coarser grid (9 million cells) and very small di¨erences in pressure were observed, while bigger di¨erences are seen in the friction coe©cient.The peak-to-peak values due to the presence of the secondary vortex are larger in the EARSM than in the kω SST model at 30 • of pitch as the predicted vortices are stronger, as seen in Fig. 24b.In this case, the main vortex is located at an azimuth of ϕ = 170 • and the secondary vortex is predicted at ϕ = 150 • .Note that an oscillation is observed very close to the secondary peak, which is due to the presence of the third vortex very close to the second one.This phenomenon is clearly captured by the EARSM.
The drag and lift predicted by both models are ¦nally compared in Table 3.As for the spheroid cases, the coe©cients were also nondimensionalized with the kω SST model.At 20 • of pitch, there is a total increase in drag of 22% when the EARSM is employed.This is due to the contributions of the pressure and friction drag components that increased by 25% and reduced by 8%, respectively.Similar  behavior is observed at 30 • of pitch where the EARSM predicted a drag 17% higher than the SST model.An increase in the lift of 12% and 9% at pitch of 20 • and 30 • , respectively, was observed.The results for C l are approximately three times higher than for the spheroid.This was expected due to the fact that the HAV has three lobes instead of one and produces more lift due to its camber.Taking this into consideration, overall, similar loadings are compared between spheroid and HAV §ows.
These results show that the turbulence model plays an important role in the prediction of this §ow.Although in the spheroid cases no large di¨erences in the loads were observed between the SST and the EARSM models, the HAV §ows appear to be more sensitive to the employed turbulence model.

Flow Topology Around a Complete Hybrid Air Vehicle Con¦guration
This subsection aims to further analyze the §ow around the complete con¦guration of a HAV.A comparison is ¦rst shown between the AL 4 and AL 5 cases of Table 1.The pressure coe©cient (C p ) at the symmetry and mid-planes are presented in Fig. 25.The solutions are practically the same, with small di¨erences close to the back of the vehicle.Grid convergence can, therefore, be assumed, with the employed URANS framework.As can also be observed in Fig. 25, favourable pressure gradients (∂p/∂s < 0) are present from the nose to 25% L, which indicates fully attached §ow.From 25% L to 80% L, a neutral pressure gradient (∂p/∂s ≈ 0) is present on both upper and lower surfaces.From 80% L, the adverse pressure gradient (∂p/∂s > 0) leads to §ow separation.The location of the transition onset on the AL body is estimated, employing empirical criteria and fully turbulent CFD results.For this, the Michel [38] and Cebeci and Smith criteria [39] are employed on streamlines extracted inside the boundary layers of the §ow solutions.These criteria are based on experimental data on §at plates and aerofoils and provide an estimate of the location of the transition onset when the local momentum thickness Reynolds number (Re θ,tr ) reaches a particular value.In the case of Michel£s criterion [38], the transition onset location takes place when Re θ,tr ≈ 2.9Re 0.4 x,tr where Re x,tr is the Reynolds number based on the distance measured from the stagnation point.In Cebeci and Smith£s correlation [39], the transition is located at the point where Re θ,tr = 1.174 1 + 22 400 Re x,tr Re 0.46 x,tr .
Figure 26 shows the momentum thickness Reynolds number (Re θ ) along a streamline of the CFD solution and the empirical transition criteria where changes in Reynolds number were considered.At Re = 3 • 10 6 , the presence of favourable pressure gradients make the boundary layer to stay laminar until approximately 80% L. When the Reynolds number is reduced, the boundary layer becomes less turbulent.The opposite e¨ect happens when the Reynolds number increases.The transition onset is moved downstream and upstream for lower  [38]; and dashdotted curves to Cebeci and Smith£s criterion [39]); and (b) transition onset and isosurfaces of reverse §ow and higher Reynolds, respectively.When the Reynolds number was reduced by a third, the transition onset was moved downstream by approximately 5% L.
When Re was increased from 3 • 10 6 to 10 • 10 6 , the transition onset moved upstream by 42% L. The results show the strong e¨ect of the Reynolds number to the nature of the boundary layer.As can be seen in the isosurfaces of the reversed §ow in Fig. 26b, when the Reynolds number is higher and, therefore, the upstream boundary layer is more turbulent, the region of the separated §ow is smaller.Conversely, at lower Reynolds number, the §ow tends to separate more.

Study of the Role of the Components on a Complete Hybrid Air Vehicle Con¦guration
To assess which components are the ones contributing the most to the overall lift and drag, Fig. 27 shows the relative changes with respect to the bare hull case (AL 1) for the 4 vehicle con¦gurations presented in Figs.2d and 3.At 10 • of pitch, as presented in Fig. 27a, the ¦ns increased the lift by more than twice and the LERX and strakes contributed to a further small increase.The drag also increased when each component was considered, but in less percentage.As Fig. 27b shows, at 20 • of pitch, the increase in loads is not as drastic as in the previous case.The ¦ns increased the lift by less than 50% and the strakes seem to provide more lift than the LERX.In this case, the penalties in drag were closer to the increase in drag than in the previous case.These results show that with the presented con¦guration of the approximated HAV, the ¦ns are the components that contribute the most to the lift and, also, drag.
Regarding the stability of the vehicle, the pitching moment (C My ) and the change of it with the pitch angle (C Mα ), calculated at the center of the volume, are presented in Fig. 28.Positive C Mα means nose up attitude and negative is  nose down.As Fig. 28b shows, both bare hull and full con¦gurations seem to be unstable as the tendency of the vehicle is to pitch up for positive changes in pitch.This should be attributed to the fact that no buoyancy force was included in the computations.Similar behavior was observed in other works in the literature that did not account for buoyancy e¨ects.This is the case of the experiments on the AKRON [40], YEZ-2A [41], and ZHIYUAN-1 [42,43] airships that also presented positive C Mα .
Nevertheless, the presence of the aerodynamic components has a stabilizing e¨ect, since there is a reduction in the pitch derivative.This was also seen by Wang et al. [42] and Freeman [40].In addition, since the ¢AL£ is not axisymmetric, at zero pitch angle, the pitching moment is not zero which di¨ers from more traditional airship con¦gurations [4144].

CONCLUDING REMARKS
The HMB2 solver was validated for the §ow around a 6:1 prolate spheroid.Within the URANS framework, the kω SST and the kω EARSM turbulence models were compared with the experimental data available in the open literature.Both models showed very good agreement in the pressure predictions.Regarding the friction coe©cients, the best agreement was obtained with the EARSM due to a lower production of eddy viscosity and better resolution of the vortices.This was observed in a comparison of the turbulent levels that showed a strong reduction of turbulence in the core of the vortices predicted by the EARSM.Analyses of the components of the stress tensor showed that the EARSM seemed to capture more spatial variations than the kω SST model.
The investigation followed by a comparison between the URANS-EARSM and SAS-SST models.The best predictions were also obtained with the URANS-EARSM, showing that to properly predict spheroid-type §ows, accounting for turbulent anisotropy is needed.Further, §ow-¦eld analyses were performed for these models.Contours of Q-criterion revealed clear di¨erences in the shape and size of the vortices, with the EARSM vortices predicted tighter and stronger.This was also observed in helicity contours.Surface friction lines showed an identical primary separation line for both URANS models, but the SAS-SST model predicted the separation at a higher azimuth angle.In addition, the EARSM predicted an earlier secondary separation line.Transition e¨ects were also investigated and it was found that the location of the transition region oscillated in time which directly a¨ected the pressure and skin friction predictions.
The e¨ect of the turbulence model on the overall loads was also presented.Di¨erences of 4% to 6% were observed in the fully turbulent cases, being the EARSM the one that predicted the highest drag.As expected, the transition model predicted much less drag (a reduction of about 30%) due to the presence of a large laminar region on the spheroid.However, for better assessment of turbulence models, experiments with integrated loads and free-transition measurements for spheroid §ows are needed.
The analysis of turbulence models ¦nished with the bare hull of a HAV.In this case, similar behavior to the spheroid£s was observed once the vehicle was pitched at 30 • .Like in spheroid §ows, the vortices were better predicted with the EARSM.The overall loads, conversely, were much more sensitive to the turbulence model in the airship case than in the spheroid.
The paper ¦nished with the analysis of the §ow around an early design of a HAV where the body vortices were identi¦ed.Due to the presence of favourable pressure gradients, the onset of transition was predicted close to the rear of the airship, at 80% L. Sensitivity studies on the Reynolds number showed a further downstream transition onset for lower Reynolds numbers and an increase of Reynolds led to an earlier transition onset.The study continued by exploring the role of each component of the airship on the aerodynamic coe©cients and the stability derivatives.The results showed that with the standard vehicle con¦guration, the ¦ns contributed the most to an increase in lift and also drag.

ACKNOWLEDGMENTS
The results were obtained using the EPSRC funded ARCHIE-West High Performance Computer (www.archie-west.ac.uk),EPSRC grant No. EP/K000586/1.The N8 HPC facility ¤Polaris¥ at Leeds University and the Chadwick cluster of the University of Liverpool were also used.Access to the systems is gratefully acknowledged as well as the ¦nancial support by the LOCATE Project.

Figure 1
Figure 1 Comparison between a spheroid and an HAV: (a) 6 : 1 prolate spheroid; and (b) hull of a representative HAV

Figure 4
Figure 4 Schematic of vortices generated around the 6 : 1 prolate spheroid at 20 • of pitch

Figure 25
Figure 25 E¨ect of the mesh density in Cp at the symmetry plane (a) and at the mid-plane (b) of the approximated HAV: 1 ¡ AL 4; and 2 ¡ AL 5

Figure 28
Figure 28 Pitching moment (CM y ) (a) and derivative (CM α ) of pitching moment (b) at the center of the volume: 1 ¡ full vehicle; and 2 ¡ bare hull

Table 1
Summary of test cases

Table 2
Drag and lift for the spheroid at 20 • of pitch , as a large laminar region was predicted with this model.Subtile di¨erences in lift are observed between fully turbulent cases.The SAS-SST model provided the lowest lift values.Conversely, the transitional model provided a much lower lift (22% lower than SAS-SST).

Table 3
Drag and lift for the HAV