A Modification to the Enhanced Correction Factor Technique for the Subsonic Wing–Body Interference Model

Abstract

In this work, the enhanced correction factor technique (ECFT) is modified for a subsonic wing–body interference model, which can consider the forces on both the lifting boxes and the body elements of an idealized airplane, termed the advanced ECFT method. A passenger aircraft model is chosen as the simulation model, and the longitudinal static aeroelasticity at the transonic situation for two-degree freedom, including the α (angle of attack) degree and ϕ (angle of horizontal tail) degree, is simulated in this paper. The corresponding CFD results are used to correct the aerodynamic influence coefficients (AIC) matrix, which is then simulated by MSC.NASTRAN. The pressure distribution results of different aircraft components received by the advanced ECFT method indicate that it is suitable for the subsonic wing–body interference model. Compared with the uncorrected linear method and the diagonal corrected method, it is generally more consistent with the CFD/CSD coupling method, not only for the lifting boxes, but also for the body elements. In addition, the aerodynamic derivative results also show good agreement with the flight test data, which solidly verifies the advance ECFT method.

1. Introduction

As early as the 1960s and 1970s, aeroelastic analysis methods based on linear aerodynamics have been widely developed [1,2]. According to the inviscid and irrotational potential flow theory, the linear aerodynamic model uses basic solutions such as vortices and doublets, to develop lifting boxes, slender elements, and interference elements, so as to obtain linearized subsonic unsteady aerodynamic data quickly. Combing this unsteady aerodynamic data with the structural flexibility and mass matrix or modal data, the aeroelastic problems can be easily and efficiently solved. This method has also been widely used in commercial software, such as NASTRAN, and is the mainstream method in aeroelastic engineering design.

The linear aerodynamic model is accurate enough for low subsonic aeroelastic problems. However, with the increase in airspeed to the transonic stage, the nonlinear effect of aerodynamics, such as shock waves, significantly changes the aerodynamic characteristics of an aircraft, which results in the linear aerodynamic model not being competent for the transonic aeroelastic analysis. In order to solve the aerodynamic nonlinearity in the aeroelastic analysis, the CFD method based on the Euler equation or N-S equation has been applied [3,4]. Combined with the structural finite element or structural modal data, a high-fidelity CFD/CSD coupled aeroelastic analysis technology has been formed. However, its huge computation consumption makes it unrealistic for massive aviation industry application.

Taking the efficiency and accuracy of the nonlinear aeroelastic analysis into account, corrections to the linear aerodynamic methods have been studied for aeroelastic analysis and design. Under a given structural modal shape, Baker et al. [5] corrected the linear DLM model using the unsteady CFD results. However, with the change in mass or stiffness, the structural modal shape is changed as well, and the unsteady CFD has to be recomputed, which makes it inefficient in engineering applications. Giesing et al. [6] used a diagonal correction matrix to correct the aerodynamic influence coefficients (AIC) matrix, which is an efficient solution for single-degree freedom correction. However, the diagonal correction matrix is not good for all modes, particularly for dissimilar modes. Hence, given the aerodynamic modes and arbitrary geometric modes, Jadic et al. [7] developed the enhanced correction factor technique (ECFT) to correct the DLM aerodynamic model, which improved the prediction of the control surface contribution. Moreno et al. [8] presented further advancement of their work, and extended the ECFT method to flutter and aeroelastic load analysis.

The ECFT method is one of the best correction methods to correct the linear DLM model in the aviation industry for the following three reasons:

• A set of experimental or CFD generated lifting forces that are used for the correction that usually already exists, which is accumulated as rigid steady aerodynamic data;
• Once determined, it can be applied repeatedly to calculations of the same condition (e.g., Mach number), irrespective of the changes in the mass and structural properties;
• A full correction matrix is conducted to correct the AIC matrix, which is suitable for muti-degree freedom correction.

However, the present ECFT method is limited to the correction of lifting boxes, which is not sufficient for the design of passenger aircraft, as the fuselage and nacelle component cannot be precisely modeled by lifting boxes, and the wing–body interaction must also be taken into consideration. The fuselage and nacelle component should be modeled as body elements, which contain two types: slender elements to simulate a body’s own motion, and interference elements to simulate the interaction with other bodies and boxes [9]. This paper presents further advancement of the present ECFT method in the field of static aeroelasticity, and makes it suitable for wing–body interference model. Meanwhile, this method is also applied to a passenger aircraft model in the transonic situation for two longitudinal degrees of freedom including α (angle of attack) degree and ϕ (angle of horizontal tail) degree. Since the AIC is corrected by the corresponding CFD results, the longitudinal static aeroelasticity simulated by MSC.NASTRAN can fit quite well with the high-fidelity CFD/CSD coupling method and the flight test.

2. Materials and Methods

2.1. Linear Aerodynamics for Wing–Body Interference Model

As for the wing–body interference DLM model, it can be simulated by singularities such as lifting boxes, slender elements, and interference elements. Downwashes can be written as follows:

$$\left\{\begin{array}{c}\underset{n_w\times 1}{\underbrace{\left[w_w\right]}}\\ \underset{\left(n_{tz}+n_{ty}\right)\times 1}{\underbrace{\left[0\right]}}\\ \underset{\left(n_{tsz}+n_{tsy}\right)\times 1}{\underbrace{\left[w_s\right]}}\end{array}\right\}=\underset{\left(n_w+n_{tz}+n_{ty}+n_{tsz}+n_{tsy}\right)\times \left(n_w+n_{tz}+n_{ty}+n_{tsz}+n_{tsy}\right)}{\underbrace{[AIC]}}\left\{\begin{array}{c}\underset{n_w\times 1}{\underbrace{\left[f_w^{\text{'}}\right]}}\\ \underset{\left(n_{tz}+n_{ty}\right)\times 1}{\underbrace{\left[{\mu }_I\right]}}\\ \underset{\left(n_{tsz}+n_{tsy}\right)\times 1}{\underbrace{\left[{\mu }_s\right]}}\end{array}\right\}$$

(1)

where:

$w_w$ = downwashes for lifting boxes;

$w_s$ = downwashes for slender elements;

$f_w^{\text{'}}$ = pressures divided by dynamic pressure along the lifting box;

${\mu }_I$ = acceleration potential for interference element doublets;

${\mu }_s$ = acceleration potential for slender element doublets;

$AIC$ = aerodynamic influence coefficients matrix;

$n_w$ = number of lifting boxes;

$n_{tz}$ or $n_{ty}$ = number of interference elements in the z (vertical) or y (lateral) orientation;

$n_{tsz}$ or $n_{tsy}$ = number of slender elements in the z (vertical) or y (lateral) orientation.

The forces related to the singularities are as follows:

$$\left\{\begin{array}{c}\underset{2n_w\times 1}{\underbrace{\left[P_w\right]}}\\ \underset{2\left(n_{tsy}+n_{tsz}\right)\times 1}{\underbrace{\left[P_s\right]}}\end{array}\right\}=\overline{q}\underset{2\left(n_w+n_{tsy}+n_{tsz}\right)\times \left(n_w+n_{tz}+n_{ty}+n_{tsz}+n_{tsy}\right)}{\underbrace{\left[SKJ\right]}}\left\{\begin{array}{c}\underset{n_w\times 1}{\underbrace{\left[f_w^{\text{'}}\right]}}\\ \underset{\left(n_{tz}+n_{ty}\right)\times 1}{\underbrace{\left[{\mu }_I\right]}}\\ \underset{\left(n_{tsz}+n_{tsy}\right)\times 1}{\underbrace{\left[{\mu }_s\right]}}\end{array}\right\}$$

(2)

and:

$$\underset{2\left(n_w+n_{tsy}+n_{tsz}\right)\times \left(n_w+n_{tz}+n_{ty}+n_{tsz}+n_{tsy}\right)}{\underbrace{\left[SKJ\right]}}=\left[\begin{array}{ccc}\underset{2n_w\times n_w}{\underbrace{\left[S_{ww}\right]}}& 0& 0\\ \underset{2\left(n_{tsy}+n_{tsz}\right)\times n_w}{\underbrace{\left[S_{sw}\right]}}& \underset{2\left(n_{tsy}+n_{tsz}\right)\times \left(n_{tz}+n_{ty}\right)}{\underbrace{\left[S_{sI}\right]}}& \underset{2\left(n_{tsy}+n_{tsz}\right)\times \left(n_{tsz}+n_{tsy}\right)}{\underbrace{\left[S_{ss}\right]}}\end{array}\right]$$

(3)

where:

$P_w$ = lifting box force and moment;

$P_s$ = slender element force and moment;

$SKJ$ = integration matrix;

$S_{ww}$ = box area for force and box area times quarter chord length for the moment;

$\left[\begin{array}{ccc}{S}_{sw}& S_{sI}& S_{ss}\end{array}\right]$ = integration matrix for slender elements, where the rows for the moment are all zero;

If only the force rows are picked up from Equation (2), it is as follows:

$$\left\{\begin{array}{c}\underset{n_w\times 1}{\underbrace{\left[F_w\right]}}\\ \underset{\left(n_{tsy}+n_{tsz}\right)\times 1}{\underbrace{\left[F_s\right]}}\end{array}\right\}=\overline{q}\underset{\left(n_w+n_{tsy}+n_{tsz}\right)\times \left(n_w+n_{tz}+n_{ty}+n_{tsz}+n_{tsy}\right)}{\underbrace{{\left[SKJ\right]}_F}}\left\{\begin{array}{c}\underset{n_w\times 1}{\underbrace{\left[f_w^{\text{'}}\right]}}\\ \underset{\left(n_{tz}+n_{ty}\right)\times 1}{\underbrace{\left[{\mu }_I\right]}}\\ \underset{\left(n_{tsz}+n_{tsy}\right)\times 1}{\underbrace{\left[{\mu }_s\right]}}\end{array}\right\}$$

(4)

And if only the moment rows are picked up from Equation (2), it is as follows:

$$\left\{\begin{array}{c}\underset{n_w\times 1}{\underbrace{\left[M_w\right]}}\\ \underset{\left(n_{tsy}+n_{tsz}\right)\times 1}{\underbrace{\left[M_s\right]}}\end{array}\right\}=\overline{q}\underset{\left(n_w+n_{tsy}+n_{tsz}\right)\times \left(n_w+n_{tz}+n_{ty}+n_{tsz}+n_{tsy}\right)}{\underbrace{{\left[SKJ\right]}_M}}\left\{\begin{array}{c}\underset{n_w\times 1}{\underbrace{\left[f_w^{\text{'}}\right]}}\\ \underset{\left(n_{tz}+n_{ty}\right)\times 1}{\underbrace{\left[{\mu }_I\right]}}\\ \underset{\left(n_{tsz}+n_{tsy}\right)\times 1}{\underbrace{\left[{\mu }_s\right]}}\end{array}\right\}$$

(5)

2.2. Advanced ECFT Method for Wing–Body Interference Model to Correct [AIC] Matrix

The first step for the advanced ECFT method is to define the downwashes $w_0$ for the linear wing–body interference model, as follows:

$$\underset{\left(n_w+n_{tz}+n_{ty}+n_{tsz}+n_{tsy}\right)\times n_m}{\underbrace{\left[w_0\right]}}=\left\{\begin{array}{c}\underset{n_w\times n_m}{\underbrace{\left[w_{0w}\right]}}\\ \underset{\left(n_{tz}+n_{ty}\right)\times n_m}{\underbrace{\left[0\right]}}\\ \underset{\left(n_{tsz}+n_{tsy}\right)\times n_m}{\underbrace{\left[w_{0s}\right]}}\end{array}\right\}$$

(6)

where $n_m$ is the number of aerodynamic modes or states.

The $f_w^{\text{'}}$ of the lifting boxes, ${\mu }_I$ of the interference elements, and ${\mu }_s$ of the slender elements due to downwashes for the linear wing–body interference model can be obtained from Equation (1):

$$\left\{\begin{array}{c}\underset{n_w\times n_m}{\underbrace{\left[f_{0w}^{\text{'}}\right]}}\\ \underset{\left(n_{tz}+n_{ty}\right)\times n_m}{\underbrace{\left[{\mu }_{0I}\right]}}\\ \underset{\left(n_{tsz}+n_{tsy}\right)\times n_m}{\underbrace{\left[{\mu }_{0s}\right]}}\end{array}\right\}=\underset{\left(n_w+n_{tz}+n_{ty}+n_{tsz}+n_{tsy}\right)\times \left(n_w+n_{tz}+n_{ty}+n_{tsz}+n_{tsy}\right)}{\underbrace{{[AIC]}^{-1}}}\cdot \underset{\left(n_w+n_{tz}+n_{ty}+n_{tsz}+n_{tsy}\right)\times n_m}{\underbrace{\left[w_0\right]}}$$

(7)

To generate aerodynamic corrections, a set of experimental or CFD generated lifting forces need to be mapped for the lifting boxes and slender elements. The aerodynamics after mapping can be stated as follows:

$$\underset{\left(n_w+n_{tsy}+n_{tsz}\right)\times n_m}{\underbrace{\left[F_1\right]}}=\left\{\begin{array}{c}\underset{n_w\times n_m}{\underbrace{\left[F_{1w}\right]}}\\ \underset{\left(n_{tsy}+n_{tsz}\right)\times n_m}{\underbrace{\left[F_{1s}\right]}}\end{array}\right\}$$

(8)

The first problem encountered in the advanced ECFT method for the wing–body interference model is that the matrix ${\left[SKJ\right]}_F$ in Equation (4) is an unfilled row, which means that after substituting Equation (8) into Equation (4), there will be more than one solution to ${\left[\begin{array}{ccc}{f}_{1w}^{\text{'}}& {\mu }_{1I}& {\mu }_{1s}\end{array}\right]}^T$. To find an appropriate solution, three assumptions are given, as follows (these three assumptions will run through the full paper and are acceptable in aviation engineering):

1.

Assumption 1, the number of interference elements is equal to the slender elements, which means$n_{tsy}=n_{ty}$ and $n_{tsz}=n_{tz}$;

2.

Assumption 2, there is one main aerodynamic mode or state in $n_m$, such as the α mode for the longitudinal problem or β (angle of sideslip angle) mode for the lateral problem;

3.

Assumption 3, the proportion of ${\mu }_{1s}/{\mu }_{1I}$ for each mode or state is equal to the ${\mu }_{0s}/{\mu }_{0I}$ for the main mode or state.

Suggest the proportion of ${\mu }_{0s}/{\mu }_{0I}$ for the main mode or state is $k_{si}$ Equation (4) can be rewritten as follows:

$$\left\{\begin{array}{c}\underset{n_w\times n_m}{\underbrace{\left[F_{1w}\right]}}\\ \underset{\left(n_{tsy}+n_{tsz}\right)\times n_m}{\underbrace{\left[F_{1s}\right]}}\end{array}\right\}=\overline{q}\left[\begin{array}{cc}\underset{n_w\times n_w}{\underbrace{{\left[S_{ww}\right]}_F}}& 0\\ \underset{\left(n_{tsy}+n_{tsz}\right)\times n_w}{\underbrace{{\left[S_{sw}\right]}_F}}& \underset{\left(n_{tsy}+n_{tsz}\right)\times \left(n_{tsz}+n_{tsy}\right)}{\underbrace{{\left[S_{sI}\right]}_F+k_{si}\left[S_{ss}\right]{}_F}}\end{array}\right]\left\{\begin{array}{c}\underset{n_w\times n_m}{\underbrace{\left[f_{1w}^{\text{'}}\right]}}\\ \underset{\left(n_{tsz}+n_{tsy}\right)\times n_m}{\underbrace{\left[{\mu }_{1I}\right]}}\end{array}\right\}=\overline{q}{\left[SKJ\right]}_F^{\text{'}}\left\{\begin{array}{c}\underset{n_w\times n_m}{\underbrace{\left[f_{1w}^{\text{'}}\right]}}\\ \underset{\left(n_{tsz}+n_{tsy}\right)\times n_m}{\underbrace{\left[{\mu }_{1I}\right]}}\end{array}\right\}$$

(9)

where, ${\left[SKJ\right]}_F^{\text{'}}$ is with filled row, and:

$$\left\{\begin{array}{c}\underset{n_w\times n_m}{\underbrace{\left[f_{1w}^{\text{'}}\right]}}\\ \underset{\left(n_{tsz}+n_{tsy}\right)\times n_m}{\underbrace{\left[{\mu }_{1I}\right]}}\end{array}\right\}=\frac{1}{\overline{q}}{\left({\left[SKJ\right]}_F^{\text{'}}\right)}^{-1}\left\{\begin{array}{c}\underset{n_w\times n_m}{\underbrace{\left[F_{1w}\right]}}\\ \underset{\left(n_{tsy}+n_{tsz}\right)\times n_m}{\underbrace{\left[F_{1s}\right]}}\end{array}\right\},\underset{\left(n_{tsz}+n_{tsy}\right)\times n_m}{\underbrace{\left[{\mu }_{1s}\right]}}=k_{si}\cdot \underset{\left(n_{tsz}+n_{tsy}\right)\times n_m}{\underbrace{\left[{\mu }_{1I}\right]}}$$

(10)

Hence:

$$\underset{\left(n_w+2n_{tsz}+2n_{tsy}\right)\times n_m}{\underbrace{\left[f_1^{\text{'}}\right]}}=\left\{\begin{array}{c}\underset{n_w\times n_m}{\underbrace{\left[f_{1w}^{\text{'}}\right]}}\\ \underset{\left(n_{tsz}+n_{tsy}\right)\times n_m}{\underbrace{\left[{\mu }_{1I}\right]}}\\ \underset{\left(n_{tsz}+n_{tsy}\right)\times n_m}{\underbrace{\left[{\mu }_{1s}\right]}}\end{array}\right\}$$

(11)

The fundamental idea behind the generation of the aerodynamic correction is to solve the system of equations as follows:

$$\underset{\left(n_w+2n_{tsz}+2n_{tsy}\right)\times \left(n_w+2n_{tsz}+2n_{tsy}\right)}{\underbrace{\left[WJJ\right]}}\cdot \underset{\left(n_w+2n_{tsz}+2n_{tsy}\right)\times \left(n_w+2n_{tsz}+2n_{tsy}\right)}{\underbrace{{[AIC]}^{-1}}}\cdot \underset{\left(n_w+2n_{tsz}+2n_{tsy}\right)\times n_m}{\underbrace{\left[w_0\right]}}=\underset{\left(n_w+2n_{tsz}+2n_{tsy}\right)\times n_m}{\underbrace{\left[f_1^{\text{'}}\right]}}$$

(12)

To make this system determinate, the null space of $\left[w_0\right]$ a matrix of rank $\left(n_w+2n_{tsz}+2n_{tsy}\right)$ is obtained as follows:

$$\underset{\left(n_w+2n_{tsz}+2n_{tsy}\right)\times \left(n_w+2n_{tsz}+2n_{tsy}\right)}{\underbrace{\left[{\mathsf{\Omega }}_0\right]}}=\left[\begin{array}{ccc}\underset{\left(n_w+2n_{tsz}+2n_{tsy}\right)\times n_m}{\underbrace{\left[w_0\right]}}& ⋮& \underset{\left(n_w+2n_{tsz}+2n_{tsy}\right)\times \left(n_w+2n_{tsz}+2n_{tsy}-n_m\right)}{\underbrace{\left[null\left(w_0\right)\right]}}\end{array}\right]$$

(13)

$\left[f_1^{\text{'}}\right]$ can also be converted to the complimentary null subspace:

$$\underset{\left(n_w+2n_{tsz}+2n_{tsy}\right)\times \left(n_w+2n_{tsz}+2n_{tsy}\right)}{\underbrace{\left[F_1^{\text{'}}\right]}}=\left[\begin{array}{ccc}\underset{\left(n_w+2n_{tsz}+2n_{tsy}\right)\times n_m}{\underbrace{\left[f_1^{\text{'}}\right]}}& ⋮& \underset{\left(n_w+2n_{tsz}+2n_{tsy}\right)\times \left(n_w+2n_{tsz}+2n_{tsy}-n_m\right)}{\underbrace{{[AIC]}^{-1}\cdot \left[null\left(w_0\right)\right]}}\end{array}\right]$$

(14)

So, the final system of Equation (12) is as follows:

$$\underset{\left(n_w+2n_{tsz}+2n_{tsy}\right)\times \left(n_w+2n_{tsz}+2n_{tsy}\right)}{\underbrace{\left[WJJ\right]}}\cdot \underset{\left(n_w+2n_{tsz}+2n_{tsy}\right)\times \left(n_w+2n_{tsz}+2n_{tsy}\right)}{\underbrace{{[AIC]}^{-1}}}\cdot \underset{\left(n_w+2n_{tsz}+2n_{tsy}\right)\times \left(n_w+2n_{tsz}+2n_{tsy}\right)}{\underbrace{\left[{\Omega }_0\right]}}=\underset{\left(n_w+2n_{tsz}+2n_{tsy}\right)\times \left(n_w+2n_{tsz}+2n_{tsy}\right)}{\underbrace{\left[F_1^{\text{'}}\right]}}$$

(15)

A unique solution for $\left[WJJ\right]$ can be obtained from Equation (15).

2.3. Advanced ECFT Method for Wing–Body Interference Model to the Correct Force and Moment

The process outlined above calculates the [WJJ] matrix, which can correct the [AIC] matrix directly. However, in some instances, the [WKK] matrix to correct the force and moment matrix might be needed, for example as an input to MSC.NASTRAN. Their difference when generating the DLM forces and moments are as follows:

$$\left\{\begin{array}{c}\underset{2n_w\times 1}{\underbrace{\left[P_w\right]}}\\ \underset{2\left(n_{tsy}+n_{tsz}\right)\times 1}{\underbrace{\left[P_s\right]}}\end{array}\right\}=\overline{q}\underset{2\left(n_w+n_{tsy}+n_{tsz}\right)\times \left(n_w+2n_{tsz}+2n_{tsy}\right)}{\underbrace{\left[SKJ\right]}}\cdot \underset{\left(n_w+2n_{tsz}+2n_{tsy}\right)\times \left(n_w+2n_{tsz}+2n_{tsy}\right)}{\underbrace{\left[WJJ\right]}}\left\{\begin{array}{c}\underset{n_w\times n_m}{\underbrace{\left[f_{0w}^{\text{'}}\right]}}\\ \underset{\left(n_{tz}+n_{ty}\right)\times n_m}{\underbrace{\left[{\mu }_{0I}\right]}}\\ \underset{\left(n_{tsz}+n_{tsy}\right)\times n_m}{\underbrace{\left[{\mu }_{0s}\right]}}\end{array}\right\}$$

(16)

$$\left\{\begin{array}{c}\underset{2n_w\times 1}{\underbrace{\left[P_w\right]}}\\ \underset{2\left(n_{tsy}+n_{tsz}\right)\times 1}{\underbrace{\left[P_s\right]}}\end{array}\right\}=\overline{q}\underset{2\left(n_w+n_{tsy}+n_{tsz}\right)\times 2\left(n_w+n_{tsy}+n_{tsz}\right)}{\underbrace{\left[WKK\right]}}\cdot \underset{2\left(n_w+n_{tsy}+n_{tsz}\right)\times \left(n_w+2n_{tsz}+2n_{tsy}\right)}{\underbrace{\left[SKJ\right]}}\left\{\begin{array}{c}\underset{n_w\times n_m}{\underbrace{\left[f_{0w}^{\text{'}}\right]}}\\ \underset{\left(n_{tz}+n_{ty}\right)\times n_m}{\underbrace{\left[{\mu }_{0I}\right]}}\\ \underset{\left(n_{tsz}+n_{tsy}\right)\times n_m}{\underbrace{\left[{\mu }_{0s}\right]}}\end{array}\right\}$$

(17)

Equations (16) and (17) should be equal to each other, which means:

$$\underset{2\left(n_w+n_{tsy}+n_{tsz}\right)\times 2\left(n_w+n_{tsy}+n_{tsz}\right)}{\underbrace{\left[WKK\right]}}\cdot \underset{2\left(n_w+n_{tsy}+n_{tsz}\right)\times \left(n_w+2n_{tsz}+2n_{tsy}\right)}{\underbrace{\left[SKJ\right]}}=\underset{2\left(n_w+n_{tsy}+n_{tsz}\right)\times \left(n_w+2n_{tsz}+2n_{tsy}\right)}{\underbrace{\left[SKJ\right]}}\cdot \underset{\left(n_w+2n_{tsz}+2n_{tsy}\right)\times \left(n_w+2n_{tsz}+2n_{tsy}\right)}{\underbrace{\left[WJJ\right]}}$$

(18)

To solve Equation (18), the [WKK] matrix and [SKJ] matrix are separated into two parts:

1

The forces part, including the forces rows (odd rows) for $\left[SKJ\right]$, recorded as ${\left[SKJ\right]}_F$, and the forces elements in odd rows and odd ranks for $\left[WKK\right]$, recorded as ${\left[WKK\right]}_F$;

2.

The moments part, including the moments rows (even rows) for $\left[SKJ\right]$, recorded as ${\left[SKJ\right]}_M$, and the moments elements in even rows and even ranks for $\left[WKK\right]$, recorded as ${\left[WKK\right]}_M$.

Equation (18) can be converted into the following:

$$\underset{\left(n_w+n_{tsy}+n_{tsz}\right)\times \left(n_w+n_{tsy}+n_{tsz}\right)}{\underbrace{{\left[WKK\right]}_F}}\cdot \underset{\left(n_w+n_{tsy}+n_{tsz}\right)\times \left(n_w+2n_{tsz}+2n_{tsy}\right)}{\underbrace{{\left[SKJ\right]}_F}}=\underset{\left(n_w+n_{tsy}+n_{tsz}\right)\times \left(n_w+2n_{tsz}+2n_{tsy}\right)}{\underbrace{{\left(\left[SKJ\right]\cdot \left[WJJ\right]\right)}_F}}$$

(19)

$$\underset{\left(n_w+n_{tsy}+n_{tsz}\right)\times \left(n_w+n_{tsy}+n_{tsz}\right)}{\underbrace{{\left[WKK\right]}_M}}\cdot \underset{\left(n_w+n_{tsy}+n_{tsz}\right)\times \left(n_w+2n_{tsz}+2n_{tsy}\right)}{\underbrace{{\left[SKJ\right]}_M}}=\underset{\left(n_w+n_{tsy}+n_{tsz}\right)\times \left(n_w+2n_{tsz}+2n_{tsy}\right)}{\underbrace{{\left(\left[SKJ\right]\cdot \left[WJJ\right]\right)}_M}}$$

(20)

In Equations (19) and (20), ${\left[SKJ\right]}_F$ and ${\left[SKJ\right]}_M$ all have unfilled rows, which means there will be no exact solution. To solve these equations, Assumption 3 is extended to Assumption 4, in order to obtain an approximate solution:

Assumption 4, the proportion of${\mu }_{0s}/{\mu }_{0I}$for each mode or state is equal to${\mu }_{0s}/{\mu }_{0I}$ for the main mode or state.

Assumption 4 can concentrate the limited error on the slender elements and interference elements, but it almost does not affect the lift boxes, which contribute the main aerodynamic force in the model.

With Assumption 4, a new ${\left[SKJ\right]}_{F_{new}}$ and ${\left(\left[SKJ\right]\cdot \left[WJJ\right]\right)}_{F_{new}}$ can be obtained, as follows:

$$\underset{\left(n_w+n_{tsy}+n_{tsz}\right)\times \left(n_w+n_{tsz}+n_{tsy}\right)}{\underbrace{{\left[SKJ\right]}_{F_{new}}}}=\left[\begin{array}{ccc}\underset{\left(n_w+n_{tsy}+n_{tsz}\right)\times n_w}{\underbrace{{\left[SKJ\right]}_F}}& ⋮& \underset{\left(n_w+n_{tsy}+n_{tsz}\right)\times \left(\left(n_w+1\right):\left(n_w+n_{tsy}+n_{tsz}\right)\right)}{\underbrace{{\left[SKJ\right]}_F}}+k_{si}\cdot \underset{\left(n_w+n_{tsy}+n_{tsz}\right)\times \left(\left(n_w+n_{tsy}+n_{tsz}+1\right):\left(n_w+2n_{tsy}+2n_{tsz}\right)\right)}{\underbrace{{\left[SKJ\right]}_F}}\end{array}\right]$$

(21)

$$\underset{\left(n_w+n_{tsy}+n_{tsz}\right)\times \left(n_w+n_{tsz}+n_{tsy}\right)}{\underbrace{{\left(\left[SKJ\right]\cdot \left[WJJ\right]\right)}_{F_{new}}}}=\left[\begin{array}{ccc}\underset{\left(n_w+n_{tsy}+n_{tsz}\right)\times \left(n_w\right)}{\underbrace{{\left(\left[SKJ\right]\cdot \left[WJJ\right]\right)}_F}}& ⋮& \underset{\left(n_w+n_{tsy}+n_{tsz}\right)\times \left(\left(n_w+1\right):\left(n_w+n_{tsy}+n_{tsz}\right)\right)}{\underbrace{{\left(\left[SKJ\right]\cdot \left[WJJ\right]\right)}_F}}+k_{si}\cdot \underset{\left(n_w+n_{tsy}+n_{tsz}\right)\times \left(\left(n_w+n_{tsy}+n_{tsz}+1\right):\left(n_w+2n_{tsy}+2n_{tsz}\right)\right)}{\underbrace{{\left(\left[SKJ\right]\cdot \left[WJJ\right]\right)}_F}}\end{array}\right]$$

(22)

However, for ${\left[SKJ\right]}_M$, the aerodynamic center of the slender element is located at the midpoint, which makes the corresponding rows zero. To modify these rows, suppose a tiny number $\epsilon$, and the aerodynamic center is $\epsilon \cdot l$ behind the midpoint, where $l$ is the length of the slender element. Hence ${\left[SKJ\right]}_{M_{new}}$ and ${\left(\left[SKJ\right]\cdot \left[WJJ\right]\right)}_{M_{new}}$ can be obtained as follows:

$$\underset{\left(n_w+n_{tsy}+n_{tsz}\right)\times \left(n_w+n_{tsz}+n_{tsy}\right)}{\underbrace{{\left[SKJ\right]}_{M_{new}}}}=\left[\begin{array}{ccc}\underset{\left(n_w+n_{tsy}+n_{tsz}\right)\times n_w}{\underbrace{{\left[SKJ\right]}_M}}& ⋮& \underset{\left(n_w+n_{tsy}+n_{tsz}\right)\times \left(\left(n_w+1\right):\left(n_w+n_{tsy}+n_{tsz}\right)\right)}{\underbrace{{\left[SKJ\right]}_M}}+k_{si}\cdot \underset{\left(n_w+n_{tsy}+n_{tsz}\right)\times \left(\left(n_w+n_{tsy}+n_{tsz}+1\right):\left(n_w+2n_{tsy}+2n_{tsz}\right)\right)}{\underbrace{{\left[SKJ\right]}_M}}\end{array}\right]$$

(23)

$$\underset{\left(n_w+n_{tsy}+n_{tsz}\right)\times \left(n_w+n_{tsz}+n_{tsy}\right)}{\underbrace{{\left(\left[SKJ\right]\cdot \left[WJJ\right]\right)}_{M_{new}}}}=\left[\begin{array}{ccc}\underset{\left(n_w+n_{tsy}+n_{tsz}\right)\times \left(n_w\right)}{\underbrace{{\left(\left[SKJ\right]\cdot \left[WJJ\right]\right)}_M}}& ⋮& \underset{\left(n_w+n_{tsy}+n_{tsz}\right)\times \left(\left(n_w+1\right):\left(n_w+n_{tsy}+n_{tsz}\right)\right)}{\underbrace{{\left(\left[SKJ\right]\cdot \left[WJJ\right]\right)}_M}}+k_{si}\cdot \underset{\left(n_w+n_{tsy}+n_{tsz}\right)\times \left(\left(n_w+n_{tsy}+n_{tsz}+1\right):\left(n_w+2n_{tsy}+2n_{tsz}\right)\right)}{\underbrace{{\left(\left[SKJ\right]\cdot \left[WJJ\right]\right)}_M}}\end{array}\right]$$

(24)

By replacing ${\left[SKJ\right]}_F$, ${\left[SKJ\right]}_M$, ${\left(\left[SKJ\right]\cdot \left[WJJ\right]\right)}_F$, and ${\left(\left[SKJ\right]\cdot \left[WJJ\right]\right)}_M$ in Equations (19) and (20) with ${\left[SKJ\right]}_{F_{new}}$, ${\left[SKJ\right]}_{M_{new}}$, ${\left(\left[SKJ\right]\cdot \left[WJJ\right]\right)}_{F_{new}}$, and ${\left(\left[SKJ\right]\cdot \left[WJJ\right]\right)}_{M_{new}}$, and ${\left[WKK\right]}_M$ can be obtained and constructed into the final $\left[WKK\right]$.

2.4. Simulation Model

A passenger aircraft model with a low horizontal tail and wing crane engine was chosen as the simulation model, as shown in Figure 1. The main parameter of this model is shown in

Table 1. Main parameters of the simulation model.

	Wing Span, m	Fuselage Length, m
Simulation model	≈35	≈39

. It has six main components, including the wing, horizontal tail, vertical tail, fuselage, nacelle, and pylon.

The wing–body interference DLM model is shown in Figure 2. The wing, horizontal tail, vertical tail, and pylon were simulated as lifting boxes. The fuselage and nacelle were simulated as body elements, including slender elements and interference elements. The width of the interference elements was unchanged from front to back, but the width of the slender elements was not. Each interference element and slender element was in both the z (vertical) and y (lateral) orientation. The number of elements is shown in

Table 2. Number of elements for the wing–body interference DLM model.

Kind of Elements	Number
Total Elements	1817
Lifting Boxes	1757
Slender Elements	30
Interference Elements	30

. MSC.NASTRAN was chosen as the simulation software.

The CFD results were used to correct the linear DLM model. Meanwhile, the CFD/CSD coupling method was used for static aeroelastic verification. The model for the CFD simulation was a cruise model, but for the CFD/CSD coupling simulation it was a jig model. The difference between them only existed in the twist angle of the wing. CFL3D was chosen as the CFD solver, and Menter’s k-Omega SST model was used. The simulation was applicable for only the longitudinal static aeroelasticity, hence half the model was used to improve the computational efficiency. The CFD mesh is shown in

Table 3. CFD grid details.

Kind of Parameters	Parameter
Total no. of points/[106]	29.2
Total no. of cells/[106]	28.9
First wall-normal layer spacing/[μm]	10
Expansion ratio	1.2
No. of wall-normal layers	40

and Figure 3.

The structural model for static aeroelastic analysis based on linear aerodynamics is shown in Figure 4. The stiffness was simulated by the beam element and the mass was simulated by the lumped mass. Auxiliary interpolation points were constructed for the interpolation of the force and displacement. Both the lumped mass points and auxiliary interpolation points are associated with beam points through rigid connections. The infinite plate spline (IPS) method was used in the interpolation between the aerodynamic model and the structural model. The inertial relief method was used for a free–free system, and the supported point was at the center of mass.

The structural model for the CFD/CSD coupling method is shown in Figure 5. A half structural model was used, corresponding to the half CFD model in Figure 3, and symmetric constraints were applied at the plane of symmetry. The beam elements and lumped mass elements were the same as the structural model in Figure 4, but the auxiliary interpolation points constructed for the interpolation of force and displacement were different when the thin plate spline (TPS) method was used. The inertial relief method [10,11] was used for the free–free system, and the supported point was at the center of mass.

Two-degree freedom, including the α degree and ϕ degree, was considered in this advanced ECFT correction, and the α degree is the main mode. The CFD results with the two-degree freedom were also simulated as a correction source to correct the linear subsonic wing–body interference DLM model.

For comparison, the following analyses were carried out:

• Rigid linear aerodynamic analysis without any correction, and static aeroelastic analysis based on linear aerodynamics without any correction;
• Rigid linear aerodynamic analysis corrected by 1 mode in the α degree, and static aeroelastic analysis corrected by 1 mode in the α degree, which is equivalent to the diagonal correction in [6];
• Static aeroelastic characteristics simulated by the high fidelity CFD/CSD coupling method, where the CFD method is based on the N-S equation.

3. Results

3.1. Rigid Pressure Difference Distribution

Three cases simulated by the CFD method were used for the advanced ECFT correction, and they had the same transonic Mach number and the same Reynold number. The detailed α and ϕ are shown in Table 4. From case 1 and case 2, an aerodynamic force for a unit α can be constructed, and from case 1 and case 3, an aerodynamic force for a unit of ϕ can be constructed. They were integrated as two modes for the advanced ECFT correction, and the [WKK] matrix was obtained in Equation (17), which was used in MSC.NASTRAN. The pressure difference distributions for a unit of α and a unit of ϕ output by different correction methods compared with the CFD results are shown in Figure 6, Figure 7 and Figure 8. For the wing component, six sections were selected along the spanwise, and for the horizontal tail, four sections were selected along the spanwise. Eta in the figures represents the relative position of the spanwise direction.

The meaning of ΔCp in figures is as follows:

$$\begin{array}{l}\mathsf{\Delta }Cp=F/\left(SQ\right)\mathrm{for}\mathrm{lifting}\mathrm{box}\\ \mathsf{\Delta }Cp=F/\left(LQ\right)\mathrm{for}\mathrm{body}\mathrm{element}\end{array}$$

(25)

where $Q$ is the dynamic pressure, $F$ is the force in each lifting box or each body element, $S$ is the area of the lifting box, and $L$ is the length of the body element. For the CFD result, $F$ is the force difference between upper and lower surface after mapping to the lifting box or body element.

Table 4. α and ϕ in the CFD simulation.

	α (°)	ϕ (°)
Case 1	0	0
Case 2	1	0
Case 3	0	1

Table 4. α and ϕ in the CFD simulation.

	α (°)	ϕ (°)
Case 1	0	0
Case 2	1	0
Case 3	0	1

As shown in Figure 6 and Figure 7, the diagonal correction was based on the α mode; hence, ΔCp received by the diagonal correction or advanced ECFT were all in good agreement with that received by CFD for a unit of α. The correction for body elements is also successful, hence the ΔCp for fuselage and nacelle after correction is also in good agreement with that received by CFD for a unit α.

However, for a unit of ϕ, as shown in Figure 8, the advanced ECFT showed obvious advantages compared with the diagonal correction. ΔCp received by the diagonal correction was close to that without any correction for a unit ϕ. For the nacelle component, ΔCp received by the diagonal correction was unreasonable. The effect of the ϕ mode on the nacelle should be particularly small. However, the ΔCp received by the advanced ECFT was in good agreement with that received by CFD, including body elements for the fuselage and nacelle.

Comparing the ΔCp received by DLM and CFD, as shown in Figure 6, Figure 7 and Figure 8, the most obvious difference was at the root of the lifting boxes, especially at the root of the horizontal tail. The DLM model simplified the aerodynamic configuration to a certain extent, and it was difficult to simulate the complex shape at the junction between the lifting boxes and the fuselage, such as the local contraction of the fuselage at the junction with the horizontal tail in this model, which led to obvious differences in the ΔCp at the root of the horizontal tail.

Figure 6. Rigid ΔCp for a unit α at Wing Component.

Figure 6. Rigid ΔCp for a unit α at Wing Component.

Figure 7. Rigid ΔCp for a unit α in the horizontal tail, fuselage, and nacelle component.

Figure 7. Rigid ΔCp for a unit α in the horizontal tail, fuselage, and nacelle component.

Figure 8. Rigid ΔCp for a unit ϕ in the horizontal tail, fuselage, and nacelle component.

Figure 8. Rigid ΔCp for a unit ϕ in the horizontal tail, fuselage, and nacelle component.

3.2. Elastic Pressure Difference Distribution

The elastic case was simulated at the same Mach number and Reynold number as the rigid case, and at a relatively large dynamic pressure, showing the relatively obvious static aeroelastic effect. For the elastic situation, with different deformations at each lifting box or each body element, the aerodynamic force would change to different levels, as shown in Figure 9, Figure 10 and Figure 11. After elastic deformation, the consistency between the diagonal correction, advanced ECFT, and CFD/CSD became poorer than the rigid situation, as shown in Figure 9 and Figure 10 for a unit of α, compared with Figure 6 and Figure 7. However, this consistency was significantly improved compared with the uncorrected results.

For a unit ϕ, as shown in Figure 11, the advanced ECFT showed obvious advantages compared with the diagonal correction. The ΔCp received by the diagonal correction was close to that without any correction. However, the ΔCp received by the advanced ECFT showed good consistency with that received by CFD/CSD.

For body elements in the fuselage and nacelle, regardless of unit α or unit ϕ, the advanced ECFT was successful, which can be seen from Figure 10 and Figure 11.

Figure 9. Elastic ΔCp for a unit α at Wing Component.

Figure 9. Elastic ΔCp for a unit α at Wing Component.

Figure 10. Elastic ΔCp for a unit α in the horizontal tail, fuselage, and nacelle component.

Figure 10. Elastic ΔCp for a unit α in the horizontal tail, fuselage, and nacelle component.

Figure 11. Elastic ΔCp for a unit ϕ in the horizontal tail, fuselage, and nacelle component.

Figure 11. Elastic ΔCp for a unit ϕ in the horizontal tail, fuselage, and nacelle component.

3.3. Elastic Aerodynamic Derivative at Different Dynamic Pressure

Elastic aerodynamic derivatives including C_Lα, C_mα, and C_mϕ are shown in Figure 12, Figure 13 and Figure 14. For C_Lα and C_mα, the results from both the diagonal correction and advanced ECFT could fit quite well with that from CFD/CSD and flight test.

Figure 12. Elastic C_Lα at different dynamic pressures.

Figure 12. Elastic C_Lα at different dynamic pressures.

Figure 13. Elastic C_mα at different dynamic pressures.

Figure 13. Elastic C_mα at different dynamic pressures.

Figure 14. Elastic C_mϕ at different dynamic pressures.

Figure 14. Elastic C_mϕ at different dynamic pressures.

However, for C_mϕ, advanced ECFT showed an obvious advantage compared with the diagonal correction. The results from the diagonal correction were close to the uncorrected results, but the results from the advanced ECFT showed good consistency with that received by the CFD/CSD and flight test.

4. Discussion

Two important functions of the aeroelastic analysis are the elastic correction of load and aerodynamic derivative. As shown in Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13 and Figure 14, the advanced ECFT achieved good results in terms of both the aerodynamic derivative and load.

Compared with the present ECFT method, the main purpose of the advanced ECFT method is to extend it from the lifting box to the subsonic wing–body interference model. As shown in Figure 6, Figure 7, Figure 8, Figure 9, Figure 10 and Figure 11, the modification is successful, where the ΔCp received by the advanced ECFT for both lifting boxes and body elements fit quiet well with that received by the CFD or CFD/CSD. This is very beneficial for the elastic correction of the load, especially for the fuselage and nacelle.

Suppose that the number of elements in the subsonic wing–body interference DLM model is n_wb, and the number of aerodynamic modes or states is n_m, there are n_wb-n_m solutions for the correct [AIC] matrix in theory. The more aerodynamic modes or states that are used, the more accurate the ECFT method, which is the advantage of the ECFT method compared with the diagonal correction method. However, when there are too many modes or states, the efficiency of the ECFT method will be greatly reduced, which is also impossible. Generally, global modes or states such as attack mode or sideslip mode, plus control surface modes or states such as horizontal tail, elevator, rudder, and aileron, are selected for the ECFT method. With the addition of the control surface modes or states in the ECFT method, the load of the control surface and the elastic control efficiency is more accurate, as shown in Figure 11 and Figure 14.

5. Conclusions

The present ECFT method is limited to the correction of lifting boxes, but is not capable for body elements including slender elements and interference elements. In the paper, the ECFT method is modified for the subsonic wing–body interference model, which is compatible with both the lifting boxes and body elements, and is termed the advanced ECFT method. Through the simulation of a transonic passenger aircraft, the pressure distribution results of different aircraft components received by the advanced ECFT method show better consistency with that received by the CFD/CSD coupling method, compared with the uncorrected linear method and the diagonal corrected method. Furthermore, the aerodynamic derivative results received by the advanced ECFT method also show good agreement with the flight test data and CFD/CSD coupling method, and have significant advantages over the uncorrected linear method and diagonal corrected method in terms of the control surface efficiency.