Verification of a low fidelity fast simulation framework through RANS simulations

Abstract

Verification and validation of simulation models are critical steps in engineering. This paper aims at verifying the suitability of reduced order aerodynamic models used in an aeroservoelastic framework designed to analyze the flight dynamics of flexible aircraft, known as the Cranfield Accelerated Aircraft Loads Model. This framework is designed for rapid assessment of aircraft configurations at the conceptual design stage. Therefore, it utilizes or relies on methods that are of relatively low fidelity for high computational speeds, such as modified strip theory coupled with Leishmann–Beddoes unsteady aerodynamic model. Hence, verification against higher order methods is required. Although low fidelity models are widely used for conceptual design and loads assessments, the open literature still lacks a comparison against higher fidelity models. This work focuses on steady-trimmed flight conditions and investigates the effect of aerodynamic wing deformation under such loads on aerodynamic performance. Key limitations of the reduced order models used, namely fuselage and interference effects, are discussed. The reasons for the overall agreement between the two approaches are also outlined.

1 Introduction

To this day, numerous simulation tools have been developed to assess aircraft non-linear flight dynamics. These include quasi-steady aerodynamic models in the form of look-up tables [4, 41] or those which utilize a complex combination of empirical and advanced models [32, 33]. For low subsonic flight, panel codes based on unsteady vortex lattice methods [31] provide a means of rapid assessment of high-altitude long-endurance aircraft. On the other hand, tools designed for high Mach number cases range from corrected doublet lattice method [39] (typically used in the loads estimations process in industry) to computationally expensive unsteady computational fluid dynamics (CFD) simulations of the whole aircraft [20]. When deployed in areas of the flight envelope where these tools are valid, designers and engineers have gained significant insight on aircraft aerodynamic performance. However, the verification and validation of these reduced order models remains a significant challenge. The authors present a comparison between the aerodynamic solver of the aeroservoelastic simulation framework [34] and a CFD results. This framework, known as Cranfield Accelerated Aircraft Loads Model (CA²LM), is designed as a MATLAB/Simulink tool for rapid simulation of numerous cases [2] as well as real-time pilot-in-the-loop simulations [28] when connected to an engineering flight simulator. As for the aerodynamics, the modeling technique relies on reduced order aerodynamic methods like modified strip theory (MST) [3] and Leishman–Beddoes unsteady models [26, 27] for key aerodynamic surfaces such as the wing, horizontal, and vertical tailplanes. On the other hand, the fuselage and nacelles are modeled using empirical ESDU methodologies [9,10,11, 13, 16].

The authors recognize that the only true method of validation for such tools is via a comparison with data obtained from a carefully tailored flight test campaign. Given the lack of access to such data and the associated costs, the discussion here focuses on a verification with the help of steady Reynolds averaged Navier–Stokes (RANS) simulations [1]. Another challenge faced when verifying such an aeroservoelastic tool is selecting a method that considers airframe flexibility. If the verification approach only considers a comparison with high fidelity CFD coupled with high fidelity computational structural mechanics [37], the computational cost of such an exercise would be too high. The aim of this research work presented in this paper is to highlight the limitations of a specific modeling approach that couples MST-based steady aerodynamic models with ESDU methods. As a result, the scope of the work discussed in this paper is limited to: (1) a comparison over a range of flight conditions and, (2) assessing the effects of wing deformation. The reader should, therefore, note that the scope of this paper is limited to the comparison of steady aerodynamic calculations; namely the RANS approach and the MST implementation. Comparisons with panel methods such as those in Refs. [18, 19, 40, 45] are considered to be beyond the scope of this paper.

This paper is structured, such that the reader is first given the details of the use case: the AX-1 large civil transport aircraft configuration. Then, a brief overview of the low fidelity method based on the classical tube and wing configuration is presented in Sect. 2. This is followed by the definition of steady CFD simulation cases in Sect. 3 that lead to a comparison with the trim solutions obtained from CA²LM in Sect. 4. In Sect. 4, the comparison is divided between undeformed and deformed airframes. Finally, Sect. 5 concludes by providing a summary of the key results.

2 Review of CA²LM aerodynamics modeling

The CA²LM framework is an aeroservoelastic simulation framework that couples structural dynamics and unsteady aerodynamics in the time domain to capture the aircraft dynamics such as structural loads. It can be used for real-time pilot-in-the-loop simulations and additional aspects such as investigating handling qualities, effect of manual controls, on flexible aircraft [29, 30] along with maneuver and gust loads [8]. The framework was developed to model a large transport aircraft known as the AX-1, as shown in Fig. 1. This framework is a MATLAB/Simulink-based tool, composed of separate blocks that carry out calculations covering the different aspects of flight mechanics. The overall architecture is shown in Fig. 2 and the reader is referred to Refs. [2, 3, 8, 28, 30, 34] for a thorough discussion of this aeroservoelastic framework.

Fig. 1

AX-1 use-case configuration

Fig. 2

CA²LM framework architecture

2.1 Modeling aerodynamic surfaces

To obtain the aerodynamic and moments steady characteristics first the wing need to be converted into a non-crank wing, this is done by employing the method presented in the ESDU 76003 [15]. With this, the aerodynamic forces and moments block in CA²LM can calculate the wing loading, that is composed of steady aerodynamic coupled with unsteady aerodynamic models.

As for the steady aerodynamic, they are implemented by the use of MST [5, 6, 42]. MST enable the calculation of the $C_\mathrm{L}$, $C_{\mathrm{D}_i}$, and $C_\mathrm{m}$ along the wingspan while assuming that the wing has no discontinuities in twist. This is done by representing a series of horseshoe vortices to calculate the downwash (w) at any point along the wingspan [3] as presented in Fig. 3 and represented by the following:

$$\begin{aligned} \begin{aligned} w(x,y)&= w_\mathrm{t}(x,y)+w_\mathrm{b}(x,y) \\&= \frac{\mathrm{d}\varGamma }{4\pi h_\mathrm{t}}(\cos \theta _1+\cos \theta _2)+\frac{\varGamma h_\mathrm{b}\mathrm{d}s}{4\pi r^3}, \end{aligned} \end{aligned}$$

(1)

Fig. 3

Representation of horseshoe vortices along the wingspan [3]

where h is the perpendicular distance of any point (x, y) to the vortex line, $\mathrm{d}\varGamma$ is the trailing vortex strength at an aerodynamic station, and $\varGamma$ is the strength of the bound vortex for some small element ($\mathrm{d}s$). By solving the downwash equation, the integrated effect of a number of trailing vortex sheets and bound vortices, while applying that the flow is tangential to the plane at the three-quarter chord line, the loading distribution may be found. As for $\varGamma$, this will depend on the local $C_\mathrm{L}$ (airfoil aerodynamics) in each aerodynamic station, and these coefficients are found from pre-programmed look-up tables (LUTs). The LUTs are generated by the use of XFOIL, a computational panel method [7] for subsonic flow, and VGK method for transonic flow [14].

On the other hand, using this simplified lifting-surface theory, the effects of compressibility may be accounted for using the Prandtl–Glauert correction factor:

$$\begin{aligned} \beta = \sqrt{1-M^2}, \end{aligned}$$

(2)

where the sweep angle ($\varLambda$) of the wing is corrected for compressibility by the following:

$$\begin{aligned} \varLambda _{\beta }=\tan ^{-1}\left( \frac{\varLambda }{\beta }\right) ; \end{aligned}$$

(3)

as for relation between the local chord at the corresponding spanwise station ($b/c_{i}$), this is also taken into account by the factor $\beta$. Regarding the variation of the local lift curve along the wingspan position, this is performed by a correction factor $K_{a_i}$, as described by DeYoung and Harper [5]. This factor may be calculated as such:

$$\begin{aligned} K_{a_i}= \frac{a_i}{2\pi /\beta }. \end{aligned}$$

(4)

Now, the local ratio $b/c_{i}$ can be corrected by multiplying it by the correction factor $1/K_{a_i}$. With this, the effects of compressibility are fully taken into account. Once that forces acting on aerodynamic surfaces are generated in the wind axis system, this needs to be transferred into the body-axis system via the application of a direction cosine matrix (DCM) with the following angles: $\theta _i$ the local twist angle, $\lambda _i$ the local sweep angle, and $\gamma _i$ the local dihedral angle. The various axes systems used in CA²LM are shown in Fig. 4.

Fig. 4

Axes systems used in the CA²LM framework

The aerodynamic model also couples the wing aerodynamics with the tail aerodynamics through a wing downwash model which interacts with the tailplane. The downwash circulation strength ($\varGamma$) is given by the following:

$$\begin{aligned} \varGamma = sV{\mathbf {A_n}}(\alpha _\mathrm{ind}-\alpha _{C_{\mathrm{L}_0}}), \end{aligned}$$

(5)

where s is the span of a wing, V is the velocity, ${\mathbf {A_n}}$ is the influence coefficients matrix for the panels used in the MST approach, and $\alpha _{C_\mathrm{L}=0}$ is zero lift angle of attack. Circulation is then evaluated through a reduced order state-space model to get the indicial angle of attack for the tailplane. This coupled implementation of the MST approach and unsteady aerodynamics modeling has been found to provide a satisfactory balance between precision and computational cost [3].

2.2 Fuselage modeling

Due to their significant contributions to forces and moments acting on the aircraft, the fuselage and engine nacelle aerodynamics are also taken into account [3]. The normal forces in the fuselage fore-body and aft-body are presented in Fig. 5 [3].

Fig. 5

Normal forces acting in the fore-body and aft-body [3]

These forces can be expressed as the summation of fore-body (fb), aft-body (ab), and carry-over lift of the fuselage in the presence of the wing (bco):

$$\begin{aligned} f_\mathrm{b}&= f_\mathrm{fb}+f_\mathrm{ab}+f_\mathrm{bco}, \end{aligned}$$

(6)

$$\begin{aligned} m_\mathrm{b}&= m_\mathrm{fb}+m_\mathrm{ab}+m_\mathrm{bco}. \end{aligned}$$

(7)

Starting with the fore-body and using empirical data and slender-body theory from the ESDU 89008 [9] and 89014 [13], it is possible to approximate the fuselage as an axisymmetric body of revolution. ESDU 89014 provides the following empirical expression for the normal force coefficient in body-axes as follows:

$$\begin{aligned} C_\mathrm{N}=C_{\mathrm{N}\alpha }\sin \alpha \cos \alpha +\frac{4}{\pi }\frac{L}{D} C_\mathrm{PL}C_\mathrm{Nc}, \end{aligned}$$

(8)

where $C_{\mathrm{N}\alpha }$ is the normal force slope obtained from ESDU 89008 [9]. The terms $C_\mathrm{PL}$ and $C_\mathrm{Nc}$ represent cross-flow effects: $C_\mathrm{PL}$ is a geometrical coefficient calculated via ESDU 77028 [17], while $C_\mathrm{Nc}$ is the cross-flow normal force coefficient and is given in the form of LUTs in ESDU 89014 as a function of Mach number and angle of attack.

To calculate moments in the fore-body, the center of pressure needs to be calculated and is given forward of the axes origin by the following:

$$\begin{aligned} X_\mathrm{cp}=X_0-\frac{(C_\mathrm{m})_0}{C_\mathrm{N}}D, \end{aligned}$$

(9)

where $X_0$ is the longitudinal location of the axes origin, D is the maximum diameter of the fuselage, and the $(C_\mathrm{m})_0$ is the pitching moment coefficient and is as follows:

$$\begin{aligned} (C_\mathrm{m})_0=(C_\mathrm{m})_{0\alpha }\sin \alpha _\mathrm{fb}\cos \alpha _\mathrm{fb}-\frac{2}{\pi } \left( \frac{L}{D}\right) ^2C_\mathrm{PL}C_\mathrm{CL}C_\mathrm{Nc}. \end{aligned}$$

(10)

The first term in Eq. (10) is the value of the pitching moment coefficient at the fore-body incident angle ($\alpha _\mathrm{fb}$) and the second term represents cross-flow effects on the pitching moment. As for $C_\mathrm{CL}$, this is given by ESDU 77028 and represents the geometric coefficients for any given cross-sectional shape.

The angle of attack of the fore-body is the sum of the elastic deformation of the fuselage and the rigid body angle of attack as presented in Fig. 6 and given by the following:

$$\begin{aligned} \alpha _\mathrm{fb}=\alpha +\alpha _\mathrm{ft}. \end{aligned}$$

(11)

Fig. 6

Angle of attack of the fore-body and aft-body of the fuselage [3]

The forces and moment contribution from the fore-body in the body-axes can now be calculated as follows:

$$\begin{aligned} f_\mathrm{fb}&= qS_\mathrm{b}\begin{bmatrix} 0 \\ 0 \\ -C_\mathrm{N} \end{bmatrix}, \end{aligned}$$

(12)

$$\begin{aligned} m_\mathrm{fb}&= qS_\mathrm{b}\begin{bmatrix} 0 \\ -C_\mathrm{N}X_\mathrm{ep} \\ 0 \end{bmatrix}, \end{aligned}$$

(13)

where $S_\mathrm{b}$ is the reference area of the fuselage given by $\pi D^2/4.$

The aft-body contributions can be estimated from slender-body theory, again approximated as an axisymmetric body and the change in normal force is given by the following:

$$\Delta C_{{\text{N}}}=\Delta C_{{{\text{N}}}\alpha }\sin \alpha \cos \alpha (1-\sin ^{0.6}\alpha )+\Delta C_{{\text{D}}_{\text{c}}}\sin ^2\alpha ,$$

(14)

$\Delta C_{\mathrm{N}\alpha }$ represents the inviscid contribution calculated according to slender-body theory and $\Delta C_{{\text{D}}_{\text{c}}}$ represents the viscous cross-flow contribution. The angle of attack for the aft-body, as presented in Fig. 6, is given by the following:

$$\begin{aligned} \alpha _\mathrm{ab}=\alpha -\alpha _\mathrm{at}. \end{aligned}$$

(15)

The center of pressure due to change in normal force is approximated in the aft-body as follows:

$$\begin{aligned} X'_\mathrm{ep}=X_0-0.5l_\mathrm{a}, \end{aligned}$$

(16)

where $l_\mathrm{a}$ is the length of the aft-body boat-tail.

Now, with this, the forces and moments about the body-axes can be calculated as follows:

$$\begin{aligned} f_\mathrm{ab}&= qS_\mathrm{b} \begin{bmatrix} 0 \\ 0 \\ -\Delta C_\mathrm{N} \end{bmatrix}, \end{aligned}$$

(17)

$$\begin{aligned} m_\mathrm{ab}&= qS_\mathrm{b}\begin{bmatrix} 0 \\ -\Delta C_\mathrm{N}X'_\mathrm{ep} \\ 0 \end{bmatrix}. \end{aligned}$$

(18)

The fuselage drag profile is estimated assuming it to be an axisymmetric body using ESDU 77028:

$$\begin{aligned} C_\mathrm{D}=\left( \frac{C_\mathrm{v}^{2/3}}{2(2\pi L/D)^{1/3}C_\mathrm{s}}\right) C_\mathrm{D}^*, \end{aligned}$$

(19)

where $C_\mathrm{D}^*$ is the body profile drag coefficient based on 2/3 of volume and the geometric coefficients $C_\mathrm{v}$ and $C_\mathrm{s}$ are calculated from the ESDU 77028 by the following equations:

$$\begin{aligned} C_\mathrm{v}&= \frac{4V}{\pi D^2L}, \end{aligned}$$

(20)

$$\begin{aligned} C_\mathrm{s}&= \frac{S}{\pi DL}. \end{aligned}$$

(21)

Now, the forces acting on the fuselage, assuming that the drag profile acts axially about the axis origin, are given as follows:

$$\begin{aligned} f_\mathrm{b}=qS_\mathrm{b}\begin{bmatrix} C_\mathrm{D} \\ 0 \\ 0 \end{bmatrix}. \end{aligned}$$

(22)

As for the nacelles, this are modeled in CA²LM as annular aerofoils using a correction for the upwash of the wing as presented in ESDU 70012 [16].

2.3 Aerodynamic interactions

The interaction between all the different parts of the aircraft is an important part of understanding the overall aerodynamic characteristics. The CA²LM framework takes into account the following interactions:

• 1. Wing–tailplane through the downwash model.

• 2. Wing–fuselage according to the ESDU 94009 [11].

Regarding the wing–fuselage interaction the lift is obtained in the outboard ($\left( C_{\ell }\right) _{\omega (f)}$) and inboard ($\left( C_{\ell }\right) _{f(\omega )}$) stations as follows:

$$\begin{aligned} (C_{\ell })_{\omega (f)}&= K_{\omega (f)}(C_{\ell _i}), \end{aligned}$$

(23)

$$\begin{aligned} (C_{\ell })_{f(\omega )}&= K_{f(\omega )}(C_{\ell _i}), \end{aligned}$$

(24)

where $K_{\omega (f)}$ and $K_{f(\omega )}$ are two factors obtained from the ESDU 91007 [12]. These effectively modify the wing lift and drag contributions calculated at various aerodynamic stations if these are within the interference distance. However, the interaction between the vertical stabilizer and the fuselage is not modeled.

3 Definition of RANS simulations

Analyzing the entire aircraft using CFD not only allows the calculation of overall aerodynamic forces and moments, but also gives an insight into the interactions between the various airframe components: like the impact of the fuselage on wing aerodynamics. Hence, the viscous effects are modeled using a Reynolds averaged Navier Stokes (RANS) solver. However, to carry out such simulations, flight conditions of interest must be clearly specified. These are defined based on the operating envelope of the AX-1 model and the conditions over which methods within the CA²LM framework are valid. Since the CA²LM framework is capable of simulating rigid and elastic airframes, CFD simulations for two configurations were carried out:

• 1. Undeformed (U) airframe over a range of flight conditions to obtain lift and drag polars (cases 1–6).

• 2. Deformed (D) airframe in flight shape at one flight condition to study the effects of airframe flexibility (case 7).

The flight conditions for performing the CFD simulations are summarized in Table 1 and Fig. 7.

Table 1 Analysis cases for underformed and deformed configurations

Fig. 7

Mission profile selected points

The value of $\alpha$ varies from $-\,3^\circ$ to $3^\circ$ for cases 1–5, and from $-\,2^\circ$ to $4^\circ$ for cases 6 and 7. The AX-1 aircraft geometry used for RANS simulations was developed based on data used in the CA²LM framework [3]. Both geometries were developed using a computer-aided design software (CATIA), and are shown in Fig. 8.

Fig. 8

Models of undeformed (cases 1–6) and deformed (case 7) cases

After defining the simulation cases and developing the aircraft models, a control volume is selected, so that a suitable mesh can be generated for CFD. The geometry of the control volume needs to be defined carefully because of the impact on the mesh generation process and the consequent results. For this study, the C-type control volume was chosen (see Fig. 9), because it was convenience for the meshing process. The size is defined in terms of the fuselage length (L). The inlet radius is equal to 5 L and aft of the aircraft the size is of 10 L. To capture the aerodynamic features affecting lift and drag in a more precise way, a refinement zone was created around the model, as shown in Fig. 10.

Fig. 9

C-type control volume of the generated mesh

Fig. 10

Density refinement zone around the aircraft model

The next step before starting the CFD analysis is to create the mesh. Therefore, a baseline mesh was defined and a grid convergence study was carried out. The parameters of the baseline mesh are given in Table 2. This base mesh contains 17.5 million elements with a mean quality value of 0.82 and a skewness value of 0.15 for only 0.003% of the elements, ensuring that the results will be accurate and that the simulation will not fail. The baseline mesh is presented in Fig. 11.

Fig. 11

Baseline mesh for AX-1 aircraft

Table 2 Baseline mesh parameters

Once the baseline mesh was defined, a grid convergence study was carried out to select the adequate meshing parameters: choosing the suitable trade-off between quality and computational cost of the simulation. The new meshes were created using the parameter in Table 2 with an increase or decrease ratio of 1.5. This study was carried out with four meshes with the number of elements ranging from around 8–26 million. And, as for the turbulence model to be implemented, several options are available. The most commonly used in the aerodynamic analysis are as follows:

• Spalart–Allmaras 1-equation model.

• $k-\epsilon$ and $k-\omega$ 2-equation model.

• SST transition 4-equation model.

The differences lay in the number of extra terms or equations added to Navier–Stokes equations. Starting with the model known as SST transition model, this is a four-equation model that combines the shear stress transport (SST) $k-\omega$ equations plus one equation for intermittency ($\gamma$) and another for the transition sources criteria ($Re_{\theta }$) [24]. As for the two-equation models, the main ones are the $k-\epsilon$ and $k-\omega$, the first one is the simplest model of turbulence in which the solution to two separate equations allows the turbulence velocity and length scales to be independently determined by the use of the transport equation for the turbulence kinetic energy (k) and its dissipation rate ($\epsilon$) [25]. The second model $k-\omega$ is different due to modifications intended to capture low Reynolds numbers effects such as compressibility and shear flow spreading by the use of the transport equation for the turbulence kinetic energy (k) and the specific dissipation rate ($\omega$) [23, 43, 44]. The Spalart–Allmaras model, on the other hand, uses only one equation to solve the kinematic eddy viscosity by solving the viscosity affected region of the boundary layer [38]. Increasing the number of equations allows the viscous effects and turbulence models to be predicted and solved more accurately, but this has the drawback of higher computational cost as presented in Fig. 12.

Fig. 12

Trade-off between computational cost and accuracy

For this study, a two-equation model, more specifically, the realizable $k-\epsilon$ with a standard wall function was selected, because it gives a suitable trade-off between computational cost and accuracy. A total of 21 simulations were run using the Cranfield University HPC system, for which each simulation consisted of 1500 iterations.

After performing all the analyses, the meshes presented a mean $y^+$ (non-dimensional unit that measures the wall distance times the shear velocity divided by the kinematic viscosity [22]) value varying between 70 and 73, indicating that the results will follow the log law [35, 36]. The actual comparison between all the meshes was done by observing changes in both the lift ($C_\mathrm{L}$) and drag ($C_\mathrm{D}$) as the mesh was further refined. Figure 13 shows how these coefficients vary with respect to computational cost.

Fig. 13

Comparison of $C_\mathrm{L}$ and $C_\mathrm{D}$ for different mesh fineness

The grid convergence study shows that as the mesh density decreases, the computational time decreases, but the precision of $C_\mathrm{L}$ and $C_\mathrm{D}$ results start to diverge. On the other hand as the mesh density increases the computational time and precision of $C_\mathrm{L}$ and $C_\mathrm{D}$ results increases accordingly. With this in mind and given that the difference in computational time between the extra fine and the fine mesh is of 136%, it was decided to keep using the parameters of the fine mesh given in Table 2 for all other meshes.

4 Comparison between CA²LM and RANS

In this section, a comparison between the aerodynamic outputs from the CA²LM framework and those obtained from the RANS simulations is presented. Moreover, this comparison is also carried out for contributions from the various airframe components. It should be noted that this study is limited to linear aerodynamics and only the lift and drag coefficients are compared. The lift coefficient ($C_\mathrm{L}$) is compared by considering the lift curve slope ($C_{\mathrm{L}_\alpha }$); this was obtained by making a linear regression due that the data are very linear. Another parameter taken into account for the $C_\mathrm{L}$ comparison is the lift coefficient at zero angle of attack ($C_{\mathrm{L}_{\alpha =0}}$). The relative differences ($\Delta _\mathrm{L}$) for both parameters are calculated as follows:

$$\begin{aligned} \Delta _\mathrm{L} = \frac{\hat{C}_{[]}-C_{[]}}{C_{[]}}, \end{aligned}$$

(25)

where $\hat{C}$ is the RANS data and C is the data from CA²LM.

Regarding the drag coefficient, starting form the traditional model:

$$C_{\text{D}} = C_{{{\text{D}}_0}}+C_{{{{\text{D}}}_{\text{i}}}},$$

(26)

where in practice $C_{\mathrm{D}_i}$ is approximately a function of $C_\mathrm{L}^2$ and conventionally the total drag coefficient is written as follows:

$$\begin{aligned} C_\mathrm{D} = C_{\mathrm{D}_0}+kC_\mathrm{L}^2, \end{aligned}$$

(27)

where k represents the lift-induced drag constant and the simple way to obtained depends on the aspect ratio (AR) and the Oswald wing efficiency factor described as follows:

$$\begin{aligned} k = \frac{1}{\pi \cdot e \cdot \mathrm{AR}}. \end{aligned}$$

(28)

Hence, the difference in drag coefficient ($\Delta _\mathrm{D}$) is done in the same way that for $\Delta _\mathrm{L}$ by comparing two values. However, having that the $C_\mathrm{D}$ depends on the parabolic model, it was decided to use instead of $C_{\mathrm{D}_0}$ that is constant. To obtain the $C_{\mathrm{D}_0}$, the values of $C_\mathrm{D}$ and $C_\mathrm{L}$ will be obtained from the RANS analysis and the data from CA²LM at 0° of the angle of attack. As for the k instead of using Eq. (28), a more exact formulation described by Howe [21] for subsonic aircraft with moderate-to-high AR is employed:

$$\begin{aligned} k = \left[ \frac{\left( 1+0.12M^6\right) }{\pi {\rm AR}}\left( 1+\frac{\left( 0.142+f(\lambda )\mathrm{AR}(10t/c)^{0.33}\right) }{\left( \cos \varLambda _{1/4} \right) ^2}+\frac{0.1\left( 3N_e+1\right) }{\left( 4+\mathrm{AR}\right) ^{0.8}}\right) \right] , \end{aligned}$$

(29)

where $N_\mathrm{e}$ is the number of engines employed, t/c is the airfoil thickness, and $f(\lambda )$ is a taper ratio function given by the following:

$$\begin{aligned} f(\lambda ) = 0.005(1+1.5(\lambda -0.6)^2). \end{aligned}$$

(30)

Taking into account that for this study the model does not have engines, that $\lambda$ is equal to 0.235 and that the airfoil employed is the NASA-SC(2)0610 with a t/c of 0.1. Table 3 presents the different value of k according to the different cases, as showed in Table 1.

Table 3 Lift-induced drag constant (k) for cases 1–5

As can be see in Table 3, the value of k is really close between all cases, because of this, an average value of 0.04349 will be employed to obtain $C_{\mathrm{D}_0}$.

In the following subsections, the comparisons for $C_\mathrm{L}$ and $C_\mathrm{D}$ for cases 1–5 are presented. These focus on the undeformed airframe and aim to verify the outputs from CA²LM over a range of flight conditions. First, overall airframe aerodynamics are considered before separate contributions from each subcomponent (wing, fuselage, HTP, and VTP) are studied. It is, therefore, possible to observe the contributions of $C_\mathrm{L}$ and $C_\mathrm{D}$ for the different components and identify areas where corrections are necessary.

4.1 Lift and drag comparison for the entire aircraft

The overall aircraft lift coefficient obtained using CA²LM compares well with the results from the RANS simulations. The largest difference is observed for case 2 where $C_{\mathrm{L}_\alpha }$ and $C_{\mathrm{L}_0}$ differ by around 8%. Figure 14a and Table 4 show that the steady aerodynamics model in CA²LM is consistently underestimating the overall lift.

With regards to the drag coefficient, the differences range between 12 and 18% where the most significant difference is observed for case 2. The overall trend is shown in Table 4 and Fig. 14b. It can be stated that the drag modeling method in CA²LM leads to significant drag overestimation and, therefore, tends to provide more conservative results.

Fig. 14

Results comparison between RANS and CA²LM

Table 4 Aircraft lift and drag coefficient comparison for cases 1–5

4.2 Fuselage lift and drag coefficient comparison

The fuselage presents the largest difference due to the limitations of the empirical methods implemented in CA²LM. Contribution to lift coefficient differs between 270% to 287% (see Table 5). In Fig. 15a, cases 2 and 4 are presented, which represent the extreme difference. Based on these trends, it can be stated that CA²LM strongly underestimates the fuselage lift compared to RANS results.

On the other hand, fuselage drag calculated by CA²LM depends entirely on the ESDU method. Table 5 shows that the difference varies between 111 and 124%, implying that CA²LM underestimates the drag generated by the fuselage by a factor of approximately 2. The causes of this discrepancy lye in the assumptions in the ESDU method, and that CA²LM does not take into account skin friction drag along the fuselage length nor the drag produced by the interaction of the fuselage and other components. This is more clearly observable in Fig. 15b where case 5 presents the greatest difference and case 1 the lowest. It should be noted that the drag difference between the cases is negligible (cases 5 and 1 differ by approximately 5%).

Fig. 15

Fuselage result comparison between RANS and CA²LM

Table 5 Fuselage lift and drag coefficient comparison for cases 1–5

4.3 Wing lift and drag coefficient comparison

Aerodynamics associated with the wing dominate the aircraft’s aerodynamics and consequently its flight dynamics. Table 6 presents the comparison between RANS and CA²LM. It shows that for all cases, the lift coefficient difference between the two solvers is below 10%, indicating that the steady aerodynamics model in CA²LM is in good agreement with the RANS results. Figure 16a shows the extreme cases, namely cases 2 and 5.

While the drag coefficient shows that at $-\,3^\circ$ of $\alpha$, the $C_\mathrm{D}$ value is practically the same for both RANS and CA²LM (Fig. 16b), the differences start to increase, primarily due to the simplified aerodynamic modeling in CA²LM. It operates under the assumption of inviscid flow, and by doing so, contributions due to viscous effects are not calculated. The results from all other cases are presented in Table 6.

Fig. 16

Wing result comparison between RANS and CA²LM

Table 6 Wing lift and drag coefficient comparison for cases 1–5

4.4 Empennage analysis

Considering the HTP contributions shown in Fig. 17a, it is evident that both CA²LM and RANS simulations produce downforce within the angle of attack range considered here, and have relatively small values of $C_{\mathrm{L}_\alpha }$. Differences between 44% and 50% are observed for the two methods (as presented in Table 7).

Now, looking into the drag coefficient comparison, Fig. 17b presents the drag polars for cases 1 and 5. It should be noted that although the HTP aerofoil section is symmetric, the drag polars for both methods are not because of the induced angle of attack from the wing downwash. Furthermore, the simplified approach to downwash modeling in CA²LM results in a significant overestimation of HTP drag. As for the differences in $C_{\mathrm{D}_0}$, Table 7 shows that these are between 32% and 40%. These differences can be attributed on how CA²LM calculates the downwash.

Fig. 17

HTP result comparison between RANS and CA²LM

Table 7 HTP lift and drag coefficient comparison for cases

Finally, the vertical tailplane (VTP) is analyzed to see its interaction with the fuselage. As stated previously, CA²LM lacks this specific interaction. It is expected that such an interaction will increase the overall drag count of the aircraft. Moreover, here, the $C_\mathrm{L}$ analysis is not considered for the VTP, because the primary force vector (lift) is perpendicular to the flow and contributes primarily through a sideforce which leads to a yawing moment. Only the $C_\mathrm{D}$ comparison is made.

Looking into the results, Fig. 18 shows that CA²LM outputs have the same values regardless of the change in angle of attack. This implies that the reduce order model does not apply any drag change with angle of attack. The difference between CA²LM and RANS varies from almost 49% to 67% and is presented in Table 8.

Fig. 18

VTP result comparison between RANS and CA²LM for $C_\mathrm{D}$

Table 8 VTP drag coefficient comparison for cases 1–5

5 Lift and drag coefficient comparison between cases 6 and 7

The comparisons made between CA²LM and RANS for cases 1–5 provide some insight into differences as flight conditions are changed. As CA²LM was designed as a tool for flight load prediction and analysis, the effects of wing deformation need to be considered. A comparison of results obtained from cases 6 (undeformed) and 7 (deformed) allows the impact of changes in wing geometry due to flight loads on aerodynamic properties to be studied. Again, the relative differences for lift and drag coefficients are used for comparing results. In this analysis, the terms $\Delta _{{\text{L}}_{\text{g}}}$ and $\Delta _{{\text{D}}_{\text{g}}}$ are introduced to evaluate the relative differences due to the geometric changes of the wing:

$$\begin{aligned} \Delta _{{\text{L}}_{\text{g}}}= & {} \frac{\dot{C}_{[]}-\ddot{C}_{[]}}{\ddot{C}_{[]}}, \end{aligned}$$

(31)

$$\Delta _{\mathrm{D}_g}= \frac{\dot{C}_{\mathrm{D}_0}-\ddot{C}_{\mathrm{D}_0}}{\ddot{C}_{\mathrm{D}_0}}.$$

(32)

Here, $\dot{C}$ represent the data of case 7 and $\ddot{C}$ the data from case 6.

5.1 Lift and drag comparison for the entire aircraft

Considering the data in Fig. 19a, it can be seen that both RANS and CA²LM give very similar results for both the deformed and undeformed cases. In both cases, the difference in $\Delta _{\mathrm{L}_\alpha }$ and $\Delta _{\mathrm{L}_0}$ are approximately 7% and 12%, respectively. This shows that both CA²LM and RANS results maintain sufficient agreement in the case of deformed wing shape. As for the difference due to the structural flexibility effects, $\Delta _{{\rm L}_{\rm g}}$ is also fairly similar. $\Delta _{\mathrm{L}_{\mathrm{g}_\alpha }}$ is approximately 1% in both methods, while $\Delta _{\mathrm{L}_{\mathrm{g}_0}}$ reaches 20% and 25% in the RANS and CA²LM results, respectively (as shown in Table 9). This indicates a slight overestimation in CA²LM when predicting changes in $C_\mathrm{L}$ due to wing flexibility.

Figure 19b and Table 9 gather the $C_\mathrm{D}$ results for the same comparison. Whilst $\Delta _\mathrm{D}$ for case 6 is 14.8%, $\Delta _\mathrm{D}$ for case 7 is smaller with a 10.9%. This suggests a greater agreement of both methods in the deformed shape. It is also shown that, once again, changes to aerodynamic coefficients are greater in CA²LM than in RANS results, with a better agreement in the case of $C_\mathrm{D}$ than the previously discussed $C_\mathrm{L}$.

Fig. 19

Comparison of overall airframe aerodynamics

Table 9 Aircraft lift ($C_{\mathrm{L}_\alpha }$ and $C_{\mathrm{L}_{\alpha =0}}$) and drag ($C_{\mathrm{D}_{\alpha =0}}$) coefficient comparison for cases 6 and 7

5.2 Wing lift and drag coefficient comparison

Looking in more detail to the performance losses due to the deformed shape of the wing, Table 10 shows that the difference between the two solvers is less than 10%. For the undeformed configuration, CA²LM overestimates the results, but in the deformed configuration, the opposite occurs: a slight underestimation in the lift produced by the wing is presented in Fig. 20a. As for the loss in lift coefficient due to the geometry change, the difference in RANS and CA²LM is more or less similar with around 30%, with a relative difference of 4% between the two methods.

The drag coefficient on the other hand does not follow the same trends as evident in Fig. 20b. At low angles of attack all cases have a relatively small difference which increases with the angle of attack. Case 6 shows a 36.9% difference, whilst case 7 is at 32.8%, as shown in Table 10. Once again, CA²LM tends to overestimate the drag when compared to RANS, but with a better convergence in the deformed configuration. Geometry changes seem to have a larger impact for CA²LM than RANS with changes at 29% and 19%, respectively.

Fig. 20

Comparison of wing aerodynamics

Table 10 Wing lift and drag coefficient comparison for cases 6 and 7

5.3 Horizontal tailplane lift and drag coefficient comparison

The HTP is analyzed to capture the impact of wing deformation on the flow hitting the tail of the aircraft (downwash). Figure 21a presents an interesting effect in which the deformed geometry in both solvers produce less downforce than the undeformed configuration, as shown in Table 11. $C_{\mathrm{L}_{\alpha =0}}$ differences for both cases between CA²LM and RANS are around 11.5%. As for $C_{\mathrm{L}_\alpha }$, there is a significant difference of 44% in case 6 and 48% in case 7 where the lift curve slope predicted by CA²LM is significantly higher for both deformed and undeformed cases. When looking into lift generated by the wing geometry change, case 7 produces less downforce than case 6, although these are at the same flight condition. This is also evident in the differences between each of the cases: RANS simulations have a difference of around 40%, whereas CA²LM shows a 53% difference.

Figure 21b and Table 11 summarize the findings for $C_\mathrm{D}$. As expected, CA²LM overestimates the drag produced when compared to RANS. Case 6 has a difference of 35.59% and case 7 around 57.67%. It should be noted that in CA²LM, the undeformed configuration produces more drag than the deformed at angles of attack below $2^\circ$. All these differences can be attributed to the change in the downwash produced by wing deformation.

Fig. 21

Comparison of HTP aerodynamics

Table 11 HPT lift ($C_{\mathrm{L}_\alpha }$ and $C_{\mathrm{L}_{\alpha =0}}$) and drag ($C_{\mathrm{D}_{\alpha =0}}$) coefficient comparison for cases 6 and 7

6 Conclusions

This study was conducted to verify the accuracy of reduced order aerodynamic models used in CA²LM—an aeroservoelastic framework used for flight loads analysis at the very early aircraft design stages. This was done by comparing the output of the various aerodynamic modules in CA²LM with results obtained from Reynolds Averaged Navier–Stokes simulations. The main motivation behind this work was to understand the suitability of the numerous reduced order aerodynamic modeling methods which have been coupled with empirical models within CA²LM. The use case for this study was the AX-1 aircraft configuration that represents a conventional large civil transport aircraft. As stated before the aim of this work is to highlight the limitations of a specific modeling approach that couples MST-based steady aerodynamic models with ESDU methods. This verification focused on two steady aerodynamic modeling aspects: (1) a comparison over a range of flight conditions and (2) assessing the effects of wing deformation. For each, the overall lift and drag coefficients are compared together with the aerodynamic contributions from individual aircraft components.

The key finding of this study was that the overall aircraft lift and drag modeling within CA²LM produces results in satisfactory agreement with the RANS simulations. However, this agreement was found to be a consequence of numerous underestimations and overestimations at the component levels. Overall drag estimation in CA²LM was found to be too conservative and mainly driven by the limitations of the MST approach used to model the various aerodynamic surfaces. Lift and drag contributions from the empirical fuselage model were found to be grossly underestimated. Comparing the results for the cases where the wing was deformed and left undeformed, the expected reduced lift coefficient due to the reduction in effective wingspan was found. However, in all cases, the primary driver behind the differences observed for the horizontal tailplane was found to be the simplified wing downwash model implemented in CA²LM. This not only affects the accuracy of aerodynamic performance, but also has a significant impact on flight dynamic modeling where delays due to wing–tailplane interactions determine critical flight loads.

To understand the root cause of the discrepancy and to compensate for these differences, a more in-depth study focusing on analysing the lift and moment distributions along the span of the lifting surfaces, the way which the CA²LM calculate the fuselage aerodynamic characteristics and the way which the downwash is calculated for the HTP are needed and form the basis of future research.