Rapid dynamic aeroelastic response analysis of the highly flexible wing with distributed propellers influence

Abstract

A rapid nonlinear aeroelastic framework for the analysis of the highly flexible wing with distributed propellers is presented, validated and applied to investigate the propeller effects on the wing dynamic response and aeroelastic stability. In the framework, nonlinear beam elements based on the co-rotational method are applied for the large-deformation wing structure, and an efficient cylinder coordinate generation method is proposed for attached propellers at different position. By taking advantage of the relatively slow dynamics of the high-aspect-ratio wing, propeller wake is modeled as a quasi-steady skewed vortex cylinder with no updating process to reduce the high computational cost. Axial and tangential induced velocities are derived and included in the unsteady vortex lattice method. For the numerical cases explored, results indicate that large deformation causes thrust to produce wing negative torsion which limits the displacement oscillation, and slipstream mainly increases the response values. In addition, an improvement of flutter boundary is found with the increase of propeller thrust while slipstream brings a destabilising effect as a result of the increment of dynamic pressure and local lift. The great potential of distributed propellers in gust alleviation and flutter suppression of such aircraft is pointed out and the method as well as conclusions in this paper can provide further guidance.

Nomenclature

A

angular acceleration

C

aerodynamic coefficient matrix

d

translational displacements

f _g

gust frequency

I ${}_\rho$

mass moment of inertia tensor

k

linear momentum

M

beam linear mass matrix

N _B

number of blades

O _TP

thrust line

p

nodal displacements

q _m

inertia loads

q _i

internal loads

q _e

external loads

(r,φ,z)

cylindrical coordinates of propeller

R _hub

propeller hub radius

R

propeller radius

${\boldsymbol{U}_\infty }$

freestream velocity

U

structural local nodal coordinate system

U _e

beam element coordinate system

w _g

gust velocity amplitude

w

angular velocity

$\alpha$

time parameter of HHT- $\alpha$ method

$\gamma$ _a

axial vortex

$\gamma$ _t

tangential vortex

Γ

circulation distribution

$\boldsymbol{\xi}$

position vector

$\boldsymbol{\rho}$

mass density

$\boldsymbol{\pi}$

angular momentum

$\boldsymbol{\theta}$

Euler angle vector of structural nodes

Ω

propeller speed

χ

skewed angle of propeller wake

1.0 Introduction

Faced with the green development requirements of high economy, low energy consumption, low noise and low pollution, the distributed electric propulsion (DEP) system using propellers is widely employed as forward power of high-altitude long-endurance (HALE) aircraft with high-aspect-ratio wings [1, 2]. Due to the highly flexible structure, the wing will experience significant structural deformation or large vibration response after being disturbed, and traditional linear analysis methods are no longer applicable. At the same time, the propeller thrust and slipstream distributed along the spanwise direction have non-negligible interference on the aerodynamic distribution, deformation distribution and dynamic response of the flexible wing, making the aeroelastic characteristics of such vehicles more complex.

In order to study the dynamic characteristics of highly flexible aircraft under large deformation, many nonlinear aeroelastic frameworks have been developed form separate research groups. An example of the nonlinear aeroelastic simulation toolboxes is UM/NAST [3–5] from University of Michigan, based on the nonlinear strain beam model and Peters two-dimensional unsteady aerodynamic model. Palacios et al. from London Imperial College developed a simulation tool called Simulation of High-Aspect-Ratio Planes (SHARP) [6–9] to investigate the dynamic characteristics of highly flexible aircraft under large deformation. Simulation of High-Aspect-Ratio Planes (SHARP) can also be used to perform nonlinear static trim, gust load calculation, stability analysis, dynamic response analysis, flutter suppression and gust mitigation. The program named ASWING [10–12] developed by Massachusetts Institute of Technology used nonlinear beams to simulate a complete conventional layout flexible aircraft. Aerodynamic calculation was conducted by a compressible modified panel model considering the deflection of control surfaces. The flight stability, flutter characteristics and dynamic response of a high-aspect-ratio flexible aircraft were studied in time domain and frequency domain respectively. Georgia Institute University developed a program called Nonlinear Aeroelastic Trim and Stability of HALE Aircraft (NATASHA) [13–15] for HALE aircraft using the nonlinear composite beam theory and the ONERA dynamic stall model. Based on this program, the effects of structural geometric nonlinearity as well as mass loads were studied on the longitudinal trim, dynamic characteristics and flight stability of a conventional layout HALE aircraft and a Helios-like full wing aircraft with distributed propellers and control surfaces, respectively.

The above investigations focus on the influence of geometric nonlinearity, and for propellers, just thrust is introduced into the framework as a follower lumped force. Concerning other propeller effects on the aeroelastic response of highly flexible aircraft, there are few related researches, and it is attracting the attention of many researchers.

Hodges [16] analysed the influence of thrust on the bending and torsion flutter of highly flexible wings without considering other effects of propellers. The propeller was modeled as a follower thrust with a certain numerical value. Results showed that the difference in flutter speed with or without propellers was up to 11%, indicating that propeller thrust had a significant impact on the flutter of highly flexible aircraft. The literature also pointed out that ignoring the slipstream and inertia effects of propellers was an important source of numerical simulation errors.

Agostinelli et al. [17] proposed a simplified and fast calculation method for the aerodynamic interference between the wing and propellers. The effect of propeller velocities along the spanwise direction was obtained by means of the blade element theory coupled with a lifting line model, and the flexible structure was modeled by linear finite beam elements. The paper indicated that propellers had a significant influence on the flow field and spanwise lift distribution of the wing.

Alba [18] applied the actuator disk theory and cylinder vortex model to investigate the aerodynamic effects of propellers on rigid and flexible wings. The method had high computational efficiency because it treated the propeller as an ideal actuator disk with an elliptic radial circulation distribution, and predetermined the propeller wake without shedding.

An [19] described a gradient-based optimisation framework for aircraft with a distributed electric propulsion system. The work used the Doublet-Lattice method to evaluate the aerodynamic forces exerted on the wing surface, an actuator disk to calculate the propeller induced velocities, and the code Toolkit for the analysis of composite structure. The framework was applied for different configurations to minimise structural weight subject to structural and aeroelastic constraints. Results showed that, despite the improvement expected for aerodynamic efficiency, larger aspect ratio wing and more electric motors led to heavier wings, which were more susceptible to flutter.

On the basis of the University of Michigan’s Nonlinear Aeroelastic Simulation Toolbox (UM/NAST), Teixeira [20–23] introduced the propeller effects including aerodynamic loads on the hub, propeller slipstream effect and inertia loads into the nonlinear aeroelastic analysis framework. Unsteady vortex lattice method was employed to calculate the wing aerodynamic force, while propellers were modeled by the lifting line theory and viscous vortex particle method. Static aerodynamic analysis showed that propeller slipstream increased wing deformation, and the dynamic response analysis indicated that the gyroscopic effect can be ignored on a typical lightweight high-aspect-ratio wing aircraft.

From the current research status, it can be seen that the influence of propellers on the aerodynamic interference, and other static responses of flexible wings has been fully studied. However, when the aircraft is disturbed during flight, further research is needed to investigate the propeller effects on the dynamic response and stability of wings undergoing large deformation. At the same time, computational cost is also an important constraint in nonlinear aeroelastic time-domain simulation considering the influence of distributed propellers.

In view of the above issues, the objectives of this paper are to establish a fast time-domain solution framework for nonlinear dynamic aeroelastic analysis of highly flexible wings with the effects of propeller thrust and slipstream, and the work focuses on the following analysis:

1. (1) Identifying the influence of geometrically nonlinear large deformation on the coupling characteristics between the wing and propellers;

2. (2) Getting insight about the effects of propeller thrust and slipstream on the dynamic response of highly flexible wings subjected to bending excitation and gust excitation;

3. (3) Investigating the importance of propeller effects on the stability of highly flexible wings under large deformation.

The remainder of this paper is organised as follows. Firstly, Section 2 provides an overview of the proposed aeroelastic analysis framework and an introduction of the models and methods in the framework. Then, in Section 3.1, a quasi-steady skewed vortex cylinder model for propeller modeling is presented and included in the unsteady vortex lattice method to consider propeller effects. Section 3.2 derives the inertial loads of co-rotational beams and gives the dynamic equilibrium equation. Section 3.3 briefly introduces the interface models including the propeller-wing structural interaction model and the spatial spline interpolation. In Section 4, the proposed aerodynamic model and the dynamic aeroelastic model are validated by comparing with literature results. In the numerical simulation of Section 5, the effects of thrust and slipstream on wing dynamic response and aeroelastic stability are analysed and compared, demonstrating again the importance of introducing propellers into the nonlinear aeroelastic analysis of highly flexible aircraft. Finally, main conclusions are drawn in Section 6.

2.0 Analysis framework

The highly flexible wing is essentially a distributed parameter system, with nonlinear changes in structural and aerodynamic response at different positions throughout the wing at each time step. In order to further investigate the influence of distributed propellers, a rapid nonlinear aeroelastic framework for the analysis of the highly flexible wing with propeller effects is established in this paper.

2.1 Overview of the framework

Figure 1 presents an overview of the established aeroelastic analysis framework. The entire analysis framework is divided into three modules, including the structural dynamics solution module, the interface module, and the aerodynamic calculation module.

Figure 1. Nonlinear dynamic aeroelastic analysis framework considering distributed propellers’ effects.

For each time step of the dynamic simulation: (1) in the structural dynamics solution module, the dynamic equation is solved and the new geometric configuration of the wing structure is determined by current values of the loads and control inputs; (2) then in the interface module, the new geometric configurations of distributed propellers are determined by current structural deformation and the aerodynamic vortex lattice grids are updated too; (3) in the aerodynamic calculation module, the aerodynamic loads of the wing and propellers are calculated and those loads are further converted to concentrated loads at structural nodes. The dynamic process continues until the predetermined simulation time is reached.

It should be noted that in order to save computational costs, the weak coupling approach between aerodynamics and structure is adopted. That is, in each time step, the aerodynamic loads of Step n+1 are determined by the structural state of Step n and is not updated during the iterative solution of the structural dynamic equation.

2.2 Introduction of models and methods in the framework

A. Aerodynamic models

For propeller aerodynamics, advanced blade element momentum theory (ABEMT) is employed to calculate the circulation distribution of propellers and skewed vortex cylinder model (SVCM) is used to obtain the propeller induced velocities. For the wing aerodynamics under propeller slipstream, unsteady vortex lattice method (UVLM) is applied with propeller induced velocities added.

Due to the propeller rotation, the inflow velocity of the blade is much larger than that of the wing. In the unsteady aerodynamic calculation, the time step required by the blade is much smaller than that required by the wing. Additionally, due to the lightweight structure of the wing, its vibration frequencies of the first few orders, such as the first-order bending frequency and the first-order twisting frequency, are much lower than the propeller rotation frequency [24]. Therefore, the propeller flow is considered quasi-steady and its wake is predetermined as a fixed vortex cylinder considering skewed inflow to save the time cost.

The circulation distribution of propellers is calculated by ABEMT, which can consider the correction of skewed inflow [25]. Different from traditional blade element momentum theory (BEMT), ABEMT not only divides blade elements along the radial direction of propeller disc, but also divides propeller disc into elements along the circumferential direction. In each element of the propeller disc, the correction of inflow velocity is considered and traditional BEMT is applied to calculate the local circulation.

B. Structural dynamic models

The high-aspect-ratio wing configuration is usually adopted in highly flexible aircraft and beam elements are often used to model the wing structure in many researches [13–16]. The spatial beams will undergo large translations and rotations but small strains during the flight, so a geometrically nonlinear structural formulation based on co-rotational (CR) beams for wing structure modeling is employed in the framework.

Within each beam element, the node large displacements are decomposed into linear elastic deformation and nonlinear rigid motion. The former is processed using linear elasticity theory in the local coordinate system, while the latter is dealt with via the co-rotational technique. The linear elastic deformation and nonlinear rigid body motion are distinguished by the co-rotation technique, which greatly simplifies the complexity of the wing geometric non-linearity. More details about the co-rotational technique can be found in the work of Crisfield [26].

To deal with the dynamic behaviour of the co-rotational beam with two nodes and six degrees of freedom per node, the inertial loads and dynamic equilibrium equations of beams are derived in the following. HHT-α algorithm is used for time integration and the algebraic non-linear equations are solved iteratively by Newton-Raphson method at each time step.

C. Interface models

Two interface models are presented in the paper including (1) the interface between structure and propeller configurations and (2) the interface between structural and aerodynamic models.

The first interface model is used to determine the new geometric configuration of distributed propellers based on current structural deformation. In order to consider the attachment of propellers, an efficient cylinder coordinate generation method is proposed for attaching propellers at different position. The expression of the propeller cylindrical coordinate system is derived based on the local nodal coordinate system of co-rotational beams. Besides, the transformation of displacements and loads between the local propeller system and the global system is completed.

The second interface model is used to convert the deformation and loads between structural and aerodynamic nodes. Because the structural and aerodynamic models use two different sets of nodes, the framework employs spatial spline interpolation to convert aerodynamic loads to concentrated loads at structural nodes, and update the aerodynamic vortex lattice nodes based on current structural deformation.

More details about each model formulation are provided in the following.

3.0 Methodology

3.1 Aerodynamic models

The framework employs ABEMT and SVCM for propeller aerodynamics and UVLM for the wing aerodynamics with propeller effects. In this section, SVCM and UVLM are presented in the following. ABEMT changes the blade element division of traditional BEMT as explained in Section 2.2 and is not presented here.

Figure 2 illustrates how these models are integrated. At a given time step, based on the new propeller configurations, the new circulation and loads of the blades are calculated by ABEMT. At the same time, based on the new propeller circulation, propeller induced velocities can be obtained from SVCM. Then, based on the new wing motion and configuration, the circulations of the wing and its wake are updated by the UVLM with propeller induced velocities considered. At last, the aerodynamic loads of the wing and propellers are output.

Figure 2. Integration of aerodynamic models.

This paper focuses on the influence of propeller slipstream on the wing aerodynamics and neglects the influence of wing on propellers installed on the leading edge. Aerodynamic interference between propellers is also ignored due to the far distance between adjacent propellers.

3.1.1 Skewed vortex cylinder model

The aerodynamic impact of propeller slipstream on the wing is mainly divided into two parts: on the one hand, the axial induced velocity component increases the local flow velocity and dynamic pressure. On the other hand, the tangential induced velocity component triggers local upwash or downwash flow. Therefore, as shown in Fig. 3, the slipstream tube model proposed by Veldhuis [27] is adopted to divide the propeller wake into two parts: axial vortex component and tangential vortex component.

Figure 3. Different components of slipstream tube mode: (a) axial vortex, (b) tangential vortex.

The induced velocity dV of the vortex element dl at any point in space can be calculated according to the Biot-Savart Law:

(1) \begin{equation}d\boldsymbol{V} = \frac{\Gamma }{{4\pi }}\frac{{\xi \times dl}}{{{{\left| \xi \right|}^3}}}\end{equation}

where Γ represents the circulation of the vortex element and $\boldsymbol{\xi}$ is the position vector. Then the axial vortex component $\gamma$ _a and the tangential vortex component $\gamma$ _t considering the skewed inflow revision [25] can be derived as follows:

(2) \begin{equation}\left\{ \begin{array}{l}{\gamma _a} = \dfrac{{{N_B}}}{{2\pi r}} \cdot \dfrac{1}{{\sqrt {1 - (1 - {{\cos }^2}\chi ){{\sin }^2}\varphi } }}\dfrac{{d\Gamma }}{{dr}}\\[16pt]{\gamma _t} = \dfrac{{\Omega {N_B}}}{{{U_\infty }}}\cos \chi \dfrac{{d\Gamma }}{{dr}}\end{array} \right.\end{equation}

where Γ(r) is the circulation distribution of propeller blades from ABEMT, N _B is the number of blades, Ω is the propeller rotation speed, ${U_\infty }$ is the velocity of freestream, and $\chi $ is the skewed angle of the wake.

In the propeller local cylindrical coordinate system (introduced in the following), the induced velocity of the tangential vortex component at point P(r _pcos $\varphi$ _p, r _psin $\varphi$ _p, z _p) can be obtained by considering Equation (2) and integrating Equation (1) along the circumference $\varphi$ , radial r and longitudinal z.

(3) \begin{equation}\boldsymbol{V}_P^t{\left( {{r_{\rm{P}}},{\varphi _{\rm{P}}},{z_P}} \right)_*} = \int\limits_{r = {R_{hub}}}^R {\frac{{{\gamma _t}}}{{4\pi }}\int\limits_0^{2\pi } {\frac{{2\left(a_*^{\prime}\sqrt c + b_*^{\prime}\sqrt a \right)}}{{\sqrt a \left(2\sqrt {ac} + b\right)}}d\varphi dr} } \end{equation}

with

(4) \begin{equation}\left\{ \begin{array}{l}\left( {a,b,c} \right) = \left( {{r^2} + r_P^2 + z_P^2 - 2r{r_P}\cos (\varphi - {\varphi _P}),2mr\cos \varphi - 2m{r_P}\cos {\varphi _P} - 2{z_P},1 + {m^2}} \right)\\[4pt]\left( {a_x^{\prime},b_x^{\prime}} \right) = (r{z_P}\cos \varphi, - r\cos \varphi )\\[4pt]\left( {a_y^{\prime},b_y^{\prime}} \right) = (r{z_P}\sin \varphi, - r\sin \varphi )\\[4pt]\left( {a_z^{\prime},b_z^{\prime}} \right) = \left( {{r^2} - r{r_P}\cos (\varphi - {\varphi _P}),rm\cos \varphi } \right)\\[4pt]m = \tan \chi \end{array} \right.\end{equation}

where the symbol ^* represents the induced velocity components of x, y and z.

Similarly, the induced velocity of the axial vortex component at point P is:

(5) \begin{equation}\boldsymbol{V}_P^a = \left( \begin{array}{l}\int\limits_{{R_{hub}}}^R {\dfrac{{{\gamma _l}}}{{4\pi }}\int\limits_0^{2\pi } {\dfrac{{ {-}\sin \varphi + \tilde r\sin {\varphi _P}}}{D}d\varphi dr} } \\[14pt]\int\limits_{{R_{hub}}}^R {\dfrac{{{\gamma _l}}}{{4\pi }}\int\limits_0^{2\pi } {\dfrac{{ - \cos \varphi + \tilde r\cos {\varphi _P} - m\tilde z}}{D}d\varphi dr} } \\[14pt]\int\limits_{{R_{hub}}}^R {\dfrac{{{\gamma _l}}}{{4\pi }}\int\limits_0^{2\pi } {\dfrac{{ - m\sin \varphi + m\tilde r\sin {\varphi _P}}}{D}d\varphi dr} } \end{array} \right)\end{equation}

with

(6) \begin{equation}\left\{\begin{array}{l} D = {D_1}\left( {{D_2} + \sqrt {1 + {m^2}} {D_1}} \right) \\[5pt] {D_1} = \sqrt {1 + {{\tilde r}^2} + {{\tilde z}^2} - 2\tilde r\cos (\varphi - {\varphi _P})} \\[5pt] {D_2} = m\cos \varphi - m\tilde r\cos {\varphi _P} - \tilde z \\[5pt] \tilde r = {r_{\textrm P}}/r,\,\,\,\,\tilde z = {z_{\textrm P}}/r \end{array} \right.\end{equation}

Since the induced velocity of the vortex components at the propeller disc is a small quantity and can be ignored, the final induced velocity of the slipstream tube model at point P is:

(7) \begin{equation}{\boldsymbol{V}_{\textrm P}} = \boldsymbol{V}_{\textrm P}^a + \boldsymbol{V}_{\textrm P}^t\end{equation}

3.1.2 Unsteady vortex lattice method with propeller effects considered

As shown in Fig. 4, the lifting surface is divided along the spanwise and chordwise directions into N elements in total. The basic solution is arranged in each vortex ring, and the collocation point is located at the quarter chord. This arrangement naturally satisfies the Kelvin condition and fluid boundary condition.

Figure 4. Lifting surface under propeller slipstream.

The grid of the lifting surface can be updated with the deformed wing structure. The structural sectional frame is considered equivalent to the body system of the aerodynamic strip to conveniently describe the motion of the collocation point, as shown in Fig. 5.

Figure 5. Motion of the non-planar vortex lattices.

Once the configuration of the deformed wing is determined, the sectional translational displacements d (d _x , d _y , d _z )^T and Euler angles $\boldsymbol{\theta}$ (θ ₁ , θ ₂ , θ ₃) corresponding to the aerodynamic strip can be obtained from structural deformation. There is a conversion relationship between the coordinates (X, Y, Z)^T in the inertial system and the corresponding coordinates (x, y, z)^T in the body system as follows:

(8) \begin{equation}\left[ {\begin{array}{*{20}{c}}x\\y\\z\end{array}} \right] = f({\theta _1},{\theta _2},{\theta _3})\left[ {\begin{array}{*{20}{c}}{X - {d_x}}\\{Y - {d_y}}\\{Z - {d_z}}\end{array}} \right]\end{equation}

with

(9) \begin{equation}f = \left[ {\begin{array}{c@{\quad}c@{\quad}c}1 {}& 0 {}& 0\\[4pt]0& {}{\cos {\theta _1}(t)} {}& {\sin {\theta _1}(t)}\\[4pt]0& {}{ {-}\sin {\theta _1}(t)} {}& {\cos {\theta _1}(t)}\end{array}} \right]\left[ {\begin{array}{c@{\quad}c@{\quad}c}{\cos {\theta _2}(t)} {}& 0 {}& { {-}\sin {\theta _2}(t)}\\[4pt]0 {}& 1 {}& 0\\[4pt]{\sin {\theta _2}(t)} {}& 0 {}& {\cos {\theta _2}(t)}\end{array}} \right]\left[ {\begin{array}{c@{\quad}c@{\quad}c}{\cos {\theta _3}(t)}& {}{\sin {\theta _3}(t)}& {}0\\[4pt]{ {-}\sin {\theta _3}(t)}& {}{\cos {\theta _3}(t)} {}& 0\\[4pt]0& {}0 {}& 1\end{array}} \right]\end{equation}

The translational velocity V ₀ of the origin of the body system in the inertia system can correspondingly be transformed to the velocity v ₀ in the body system:

(10) \begin{equation}{\boldsymbol{v}_0} = f({\theta _1},{\theta _2},{\theta _3}){\boldsymbol{V}_0}\end{equation}

Then, in the body system, the velocity of the collocation point j is equal to:

(11) \begin{equation}\boldsymbol{v}_{_j}^{motion} = - ({\boldsymbol{v}_0} + {\textbf{w}_{\boldsymbol{b}}} \times {\boldsymbol{\xi}_j})\end{equation}

where $\boldsymbol{\xi}$ _j is the position vector of the collocation point and ${\textbf{w}_b}$ is the angular velocity of the body system relative to the inertial system.

The trailing vortex rings are shed from the trailing edge, and the circulation of the new is equal to that of the trailing vortex ring at the previous moment, namely:

(12) \begin{equation}\Gamma _W^t = \Gamma _{T.E.}^{t - \Delta t}\end{equation}

Then, the induced velocity ${\boldsymbol{v}}_{\rm{j}}^{{\rm{wake}}}$ of the wake at the collocation point j can be calculated using the Biot-Savart Law based on Equation (12).

To satisfy the non-penetration boundary condition, the induced velocity component in the normal direction n _j at the collocation point j is zero. The equations of N vortex rings can be written in the matrix form as:

(13) \begin{equation}\left[ {\begin{array}{*{20}{c}}{{\Gamma _1}}\\{{\Gamma _2}}\\ \vdots \\{{\Gamma _N}}\end{array}} \right] = {\textbf{C}^{ - 1}}\left[ {\begin{array}{*{20}{c}}{\boldsymbol{v}_1^T{\boldsymbol{n}_1}}\\{\boldsymbol{v}_2^T{\boldsymbol{n}_2}}\\ \vdots \\{\boldsymbol{v}_N^T{\boldsymbol{n}_N}}\end{array}} \right]\end{equation}

with

(14) \begin{equation}{\boldsymbol{v}_j} = {\boldsymbol{U}_\infty } + \boldsymbol{v}_j^{wake} + \boldsymbol{v}_j^{motion} + \boldsymbol{v}_j^{prop}\quad j = 1,2, \cdots, N\end{equation}

where C is the aerodynamic coefficient matrix determined by the bound vortex rings and v _j is total the induced velocity at the collocation point j, which includes the freestream velocity U _∞, the motion velocity ${\boldsymbol{v}}_{\rm{j}}^{{\rm{motion}}}$ from Equation (11), the induced velocity ${\boldsymbol{v}}_{\rm{j}}^{{\rm{wake}}}$ of wing wake from Equation (12), and the induced velocity ${\boldsymbol{v}}_{\rm{j}}^{{\rm{prop}}}$ of slipstream from Equation (7).

Due to the predetermined propeller wake, the induced velocity will be distorted when the vortex lattice is near the slipstream boundary. To solve the problem, the average induced velocity at the four corners of the vortex lattice is used to replace the induced velocity at the collocation point, which can effectively reduce the result error.

3.2 Structural dynamic models

In this section, the highly flexible wing structure is modeled by nonlinear co-rotational beams, and the inertial loads of beams and the dynamic equilibrium equation are introduced. The beams are assumed to be homogeneous, isotropic and linearly elastic. The beams are assumed to be homogeneous, isotropic and linearly elastic. Nodes are located in the line of centroids of the cross-section, which is equivalent to the elastic axis under the assumption.

As shown in Fig. 6, the beam deformed configuration is described through the local nodal coordinate systems. The left and right nodal coordinate systems of the beam are ${\boldsymbol{U}}_{\rm{1}}$ and ${{\boldsymbol{U}}_{\rm{2}}}$ , respectively.

Figure 6. Local nodal coordinate systems of CR beams.

At a certain moment, ${{\boldsymbol{U}}_{\rm{1}}}$ and ${{\boldsymbol{U}}_{\rm{2}}}$ in the deformed configuration can be determined by rotating the initial ${{\boldsymbol{U}}_{{\rm{1,0}}}}$ and ${{\boldsymbol{U}}_{{\rm{2,0}}}}$ through the Euler angle vector $\boldsymbol{\theta}$ _i of the corresponding structural nodes:

(15) \begin{equation} {{{\boldsymbol{U}}_i} = f({\boldsymbol\theta _i}){{\boldsymbol{U}}_{i,0}}}\quad {}{i = 1,2}\end{equation}

where f is the same as Equation (9).

The linear momentum k and the angular momentum $\boldsymbol{\pi}$ of beams per unit length can be defined by the following equations:

(16) \begin{equation}\left\{ \begin{array}{l}\boldsymbol{k} = A\rho \dot{\boldsymbol{d}}\\[4pt] \boldsymbol{\pi} = \boldsymbol{U}{\boldsymbol{I}_\rho }{\boldsymbol{U}^T}\textbf{w}\end{array} \right.\end{equation}

where A is the cross-section aera of the beam, $\rho$ is the mass density, $\dot{\boldsymbol{d}}$ is the translational velocity of the centre of mass of the section, U is the local structural frame of the section, w is the angular velocity of the section and I ${}_\rho$ is the mass moment of inertia tensor, defined as:

(17) \begin{equation}{\boldsymbol{I}_\rho } = \rho \left[ {\begin{array}{c@{\quad}c@{\quad}c}{{I_1} + {I_2}} & {}0 {} & 0\\[4pt]0 {} &{{I_1}} {} &0\\[4pt]0 {} &0 {} & {{I_2}}\end{array}} \right]\end{equation}

where I ₁ and I ₂ are the principal moments of inertia of the cross-section.

The inertia loads of the cross-section ${\boldsymbol{q}_{section}}$ can be obtained by taking the time derivative of Equation (16):

(18) \begin{equation}{\boldsymbol{q}_{section}} = \frac{d}{{dt}}\left( {\begin{array}{c}\boldsymbol{k}\\[4pt]\boldsymbol{\pi }\end{array}} \right) = \left( {\begin{array}{c}{A\rho \ddot{\boldsymbol{d}}}\\[4pt]{\boldsymbol{U}{\boldsymbol{I}_\rho }{\boldsymbol{U}^{\textrm{T}}}\textbf{A} + \textbf{w} \times \boldsymbol{U}{\boldsymbol{I}_\rho }{\boldsymbol{U}^{\textrm{T}}}\textbf{w}}\end{array}} \right)\end{equation}

where A is the angular acceleration of the section.

For the beam with two nodes, the shape functions N ₁(x) = (l-x)/l and N ₂(x) = x/l are used to interpolate sectional displacements through nodal displacements. Combined with the principle of virtual work, the equivalent nodal inertial loads ${\boldsymbol{q}_m}$ of the beam can be obtained as:

(19) \begin{equation}{\boldsymbol{q}_m} = \int_0^l {\left( {\begin{array}{c@{\quad}c}{{N_1}(x)\textbf{I}} {}& \textbf{0}\\[4pt]\textbf{0}& {}{{N_1}(x)\textbf{I}}\\[4pt]{{N_2}(x)\textbf{I}}& {}\textbf{0}\\[4pt]\textbf{0} {}& {{N_2}(x)\textbf{I}}\end{array}} \right)} \times \left( {\begin{array}{*{20}{c}}{A\rho \ddot{\boldsymbol{d}}}\\[4pt]{\boldsymbol{U}{\boldsymbol{I}_\rho }{\boldsymbol{U}^{\textrm{T}}}\textbf{A} + \textbf{w} \times \boldsymbol{U}{\boldsymbol{I}_\rho }{\boldsymbol{U}^{\textrm{T}}}\textbf{w}}\end{array}} \right)dx\end{equation}

where l is the length of the beam, I is the identity matrix, 0 is the zero matrix.

Reference (26) provided the specific expression of the inertial loads q _m of the nonlinear beam as:

(20) \begin{equation}{\boldsymbol{q}_m} = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}\textbf{I} {}&&&\\& {}{{\boldsymbol{U}_e}} {}& &\\& {}& {}\textbf{I} {}&\\&&& {}{{\boldsymbol{U}_e}}\end{array}} \right]\left(\textbf{M}{\ddot{\boldsymbol{p}}}+ \frac{l}{{12}}\left[ {\begin{array}{*{20}{c}}{{\textbf{0}_{3 \times 1}}}\\[4pt]{3\boldsymbol{S}({\textbf{w}_1}){\boldsymbol{I}_\rho }{\textbf{w}_1} + \boldsymbol{S}({\textbf{w}_1}){\boldsymbol{I}_\rho }{\textbf{w}_2} + \boldsymbol{S}({\textbf{w}_2}){\boldsymbol{I}_\rho }{\textbf{w}_1} + \boldsymbol{S}({\textbf{w}_2}){\boldsymbol{I}_\rho }{\textbf{w}_2}}\\[4pt]{{\textbf{0}_{3 \times 1}}}\\[4pt]{\boldsymbol{S}({\textbf{w}_1}){\boldsymbol{I}_\rho }{\textbf{w}_1} + \boldsymbol{S}({\textbf{w}_1}){\boldsymbol{I}_\rho }{\textbf{w}_2} + \boldsymbol{S}({\textbf{w}_2}){\boldsymbol{I}_\rho }{\textbf{w}_1} + 3\boldsymbol{S}({\textbf{w}_2}){\boldsymbol{I}_\rho }{\textbf{w}_2}}\end{array}} \right]\right)\end{equation}

where M is the linear mass matrix of the beam element, U _e is the local element frame of beams and it can be found by averaging the left and right nodal coordinate systems U ₁ and U ₂, respectively. w _i is the angular vector of node i defined in U _e , $\boldsymbol{p}$ is nodal displacements, which include translational components and rotational components, S (·) is an algebraic operation matrix given by:

(21) \begin{equation}\boldsymbol{S}(\textbf{w}) = \left[ {\begin{array}{c@{\quad}c@{\quad}c}0 & {}{ - {\textrm{w}_z}} & {}{{\textrm{w}_y}}\\[4pt]{{\textrm{w}_z}} & {}0 & {}{ - {\textrm{w}_x}}\\[4pt]{ - {\textrm{w}_y}} & {}{{\textrm{w}_x}} & {}0\end{array}} \right]\end{equation}

The dynamic equilibrium is imposed according to the following equation by HHT- $\alpha$ method, which is one of the most commonly used time integration schemes:

(22) \begin{equation}(1 + \alpha )({\boldsymbol{q}_{i,n + 1}} - {\boldsymbol{q}_{e,n + 1}}) - \alpha ({\boldsymbol{q}_{i,n}} - {\boldsymbol{q}_{e,n}}) + {\boldsymbol{q}_{m,n + 1}} = 0\end{equation}

Where the inertia terms q _m are related to the end-point t _n+1 from Equation (20), the internal force terms q _i and the external force terms q _e are related to some $\alpha$ -point between t _n and t _n+1.

Newton-Raphson method is used to solve such a system of algebraic nonlinear equations. The dynamic tangent stiffness matrix can be obtained by computing the variation of the internal forces q _i,n+1 and inertial forces q _m,n+1. More details about the dynamics of CR beams can be found in Ref. [26].

3.3 Interface models

In this section, the propeller-wing structural interaction model and the spatial spline interpolation are presented. The former is used for the interface between structure and propeller configurations and the latter is used for the interface between structural and aerodynamic models.

3.3.1 Propeller-wing structural interaction model

The position and direction of distributed propellers change in real time with the wing deformation. In order to convert displacements and loads between the wing and propellers, the expression of the local cylindrical coordinate systems of propellers is derived to complete the unification of the inertial system, the structural nodal coordinate system, and the propeller cylindrical coordinate system.

In the inertia system, as shown in Fig. 7, for a propeller propulsion system typically installed on the leading edge, the origin of the propeller cylindrical coordinate system is located at the centre point O _p, and the polar coordinate r _p(r _p, $\varphi$ _p) is located in the rotation plane. The direction of z _p is perpendicular to the rotation plane and opposite to the thrust line O _TP. The structural nodal coordinate system at the installation section is U _i (U _ix , U _iy, U _iz ), and the coordinate of the corresponding node is O _i .

Figure 7. Propeller cylinder coordinate system.

In the local nodal coordinate system, the coordinate of the propeller centre point is ${\boldsymbol{O}}_{\rm{p}}^{{\rm{\;l}}}$ , and the direction of thrust line is ${\boldsymbol{O}}_{{\rm{TP}}}^{{\rm{\;l}}}$ (the superscript l representing the quantity in the local system). During the structural motion process, the local quantity ${\boldsymbol{O}}_{\rm{p}}^{{\rm{\;l}}}$ and ${\boldsymbol{O}}_{{\rm{TP}}}^{{\rm{\;l}}}$ can be converted to the global quantity O _p and O _TP by following equations:

(23) \begin{equation}{\boldsymbol{O}_{\textrm{p}}} = {\boldsymbol{O}_i} + {\boldsymbol{U}_i} \cdot \boldsymbol{O}_{\textrm{p}}^{\textrm{l}}\end{equation}

(24) \begin{equation}{\boldsymbol{O}_{\textrm{TP}}} = {\boldsymbol{U}_i} \cdot \boldsymbol{O}_{\textrm{TP}}^{\textrm{l}}\end{equation}

Similarly, for any point H(x _H, y _H , z _H) in the inertial system, its coordinates ( ${\boldsymbol{r}}_{\rm{H}}^{{\rm{\;l}}}$ , ${\boldsymbol{Z}}_{\rm{H}}^{{\rm{\;l}}}$ ) in the propeller cylindrical coordinate system can be obtained from following equations:

(25) \begin{equation}\boldsymbol{z}_{\rm{H}}^{\textrm{l}} = - {(\textbf{H} - {\boldsymbol{O}_{\textrm{p}}})^\textrm{T}} \cdot {\boldsymbol{O}_{\textrm{TP}}}\end{equation}

(26) \begin{equation}\boldsymbol{r}_{\rm{H}}^{\textrm{l}} = - (\textbf{H} - {\boldsymbol{O}_{\textrm{p}}}) \times {\boldsymbol{O}_{\textrm{TP}}}\end{equation}

Reference (22) pointed out that for the investigated propellers with relatively low rotating speed and small mass, gyroscopic effects had little influence on the aeroelastic response and stability of high-altitude long-endurance aircraft. Therefore, the gyroscopic effects of the propulsion system are ignored. The thrust T _p and torque M _p of propellers can be obtained by ABEMT, and they can be converted to the structural nodes as F _i and M _i in the inertia system as follows:

(27) \begin{equation}{\boldsymbol{F}_i} = {T_{\textrm{p}}} \cdot {\boldsymbol{O}_{\textrm{TP}}}\end{equation}

(28) \begin{equation}{\boldsymbol{M}_i} = - {T_{\textrm{p}}}({\boldsymbol{O}_i} - {\boldsymbol{O}_{\textrm{p}}}) \times {\boldsymbol{O}_{\textrm{TP}}} + {M_{\textrm{p}}}{\boldsymbol{O}_{\textrm{TP}}}\end{equation}

3.3.2 Spatial spline interpolation

Spatial spline interpolation is used to complete the displacement and load conversion between structural nodes and aerodynamic nodes. For structural nodes, the translational and rotational displacements are d _s and $\boldsymbol{\theta}$ _s respectively. Assuming the rigid behaviour in the chordwise direction, the translational displacement d _a of aerodynamic nodes can be obtained by:

(29) \begin{equation}{\boldsymbol{d}_{\textrm{a}}} = {\boldsymbol{d}_{\textrm{s}}} + f({\boldsymbol\theta _s}){\boldsymbol\xi _{\textrm{a}}} = \boldsymbol{G}{\boldsymbol{d}_s}\end{equation}

where $\boldsymbol{\xi}$ _a is the vector from the aerodynamic node to the structural node of corresponding section, f is the same as Equation (9), and G is the displacement interpolation matrix, of which the specific expression is shown in Ref. [28].

Let the structural virtual displacement be $\delta$ d _s, the aerodynamic virtual displacement be $\delta$ d _a, the aerodynamic load be F _a and the structural load be F _s. According to the virtual work principle, there is:

(30) \begin{equation}\delta \boldsymbol{d}_{\textrm{a}}^\textrm{T}{\boldsymbol{F}_{\textrm{a}}} = \delta \boldsymbol{d}_{\textrm{s}}^\textrm{T}{\boldsymbol{F}_{\textrm{s}}}\end{equation}

According to Equation (29), the displacement increment relationship between structural and aerodynamic nodes can be obtained:

(31) \begin{equation}\delta {\boldsymbol{d}_{\textrm{a}}} = \boldsymbol{G}\delta {\boldsymbol{d}_{\textrm{s}}}\end{equation}

By bring Equation (31) into Equation (30), we get:

(32) \begin{equation}\delta \boldsymbol{d}_{\textrm{s}}^\textrm{T}{\boldsymbol{G}^\textrm{T}}{\boldsymbol{F}_{\textrm{a}}} = \delta \boldsymbol{d}_{\textrm{s}}^\textrm{T}{\boldsymbol{F}_{\textrm{s}}}\end{equation}

According to the arbitrariness of virtual displacement, the interpolation relationship between aerodynamic loads and structural loads can be obtained as follows:

(33) \begin{equation}{{\bf{F}}_{\textrm{s}}} = {\boldsymbol{G}^\textrm{T}}{\boldsymbol{F}_{\textrm{a}}}\end{equation}

4.0 Validation studies

The model validation is carried out in this section for the proposed aeroelastic analysis framework of the highly flexible wing considering propeller effects. On the one hand, aerodynamic models are validated for two configurations: the isolated propeller and the rigid wing coupled with propellers. On the other hand, for the highly flexible wing attached with distributed propellers, the static aeroelastic deformation and the dynamic aeroelastic response are validated respectively.

4.1 Validation of the isolated propeller

As a first check for the implementation of propeller aerodynamics, the comparison of thrust with results from Ref. [29] is performed. The propeller information and calculation conditions are summarised in Table 1. The blade is divided into 17 elements along 0.2R∼R interval 0.05R, and the propeller speed is set to $2,000 \sim 7,000$ rpm.

Table 1. Parameters of the isolated propeller

Figure 8 presents the comparison curves between thrust and propeller speed, which shows the thrust increases nonlinearly with propeller speed. It can be seen that results of this paper are in good agreement with the numerical and experimental results in Ref. [29], and the relative error of results is within 1.5%.

Figure 8. Comparison of thrust with propeller speed.

The effect of propeller on the flow field is also compared. The propeller speed is set to 4,600 rpm, and the calculation conditions such as flow velocity and density are unchanged. Figure 9 gives the distribution of the total axial velocity and the tangential induced velocity at one radius behind the propeller. It can be seen that the total axial velocity is symmetrically distributed along the rotation centre, and the tangential induced velocity is approximately antisymmetric. As discussed in Ref. [27], the axial velocity increases the local dynamic pressure, while the tangential velocity changes the local angle-of-attack with positive effect at the side where blades go up and negative at the other side. It can also be observed that results in this paper are generally consistent with the CFD time-averaged results in Ref. [30], proving that the model proposed can effectively obtain the induced effect of propeller slipstream.

Figure 9. Distribution of the (a) total axial velocity and (b) tangential induced velocity.

4.2 Validation of the propeller-wing unsteady aerodynamics

In the aerodynamic model proposed, the propeller wake near the wing is assumed to be quasi-steady with no updating process to save the computing cost. In order to validate the rationality of the hypothesis, a propeller-wing system with simple geometry is analysed and the lift coefficient is compared with that calculated by UVLM method for both propeller and wing from Ref. [24]. The parameters of the system are given in Table 2 and the wing is rigid with five propellers attached at 1m intervals along the spanwise direction and located 0.2m ahead of and 0m upon the leading edge, as shown in Fig. 10. To focus on the unsteady effect, the wing is forced to undergo the heaving motion at a low frequency of 0.6Hz and an amplitude of 0.3m.

Table 2. Parameters of the propeller and wing in aerodynamic study

Figure 10. Schematic diagram of the propeller-wing system in aerodynamic study.

According to the convergence tests in Ref. [24], the number of spanwise panels is 48 and the number of chordwise panels is 14. If the former is too small, the collocation points on the wing will be outside the propeller slipstream, which will lead to a decrease in lift. At the same time, as described in Section 3.1.2, a processing method is applied to reduce the calculation error by using the averaged value of the induced velocity at four corners of the vortex lattice instead of the induced velocity at the collocation point. Additionally, to express the high-speed-rotating wake-updating of the propeller, a small dt = 0.002s is employed for the comparison.

Figure 11 shows the time history of the wing lift coefficient in 2 seconds. Results of the clean wing are calculated by UMLV proposed in Section 3.1.2 without considering the propeller induced velocities. Due to the presence of propellers, the oscillation amplitude of the wing lift coefficient has increased significantly compared to the clean wing. Compared with literature results, the proposed method has a certain accuracy in the unsteady calculation of low dynamic frequency, and can effectively obtain the aerodynamic interference effect of propellers on the highly flexible wing.

Figure 11. Comparison of the wing lift coefficient in unsteady heaving motion at frequency = 0.6Hz and amplitude = 0.3m.

4.3 Validation of the aeroelastic response

This section mainly validates the characteristics of aeroelasticity considering the influence of propeller slipstream from two aspects: static deformation and dynamic response. The high-aspect-ratio wing [31, 32] is selected as the validation object with ten propellers symmetrically distributed at 3/8, 4/8, 5/8, 6/8 and 7/8 of the half-span wing as shown in Fig. 12. All propellers have a chordwise installation distance of 0.5m from the leading edge, a vertical installation distance of 0m, an installation angle of 180°, and a positive direction of rotation (the outside of the wing is the upwash area). The freestream velocity is 25m/s. The air density is 0.0889kg/m³. Other parameters of the propeller and wing are shown in Table 3.

Figure 12. Schematic diagram of the propeller-wing system in aeroelastic study (half model).

Table 3. Parameters of the propeller and wing in aeroelastic study

Based on the convergence study, the structural model is divided into 32 nonlinear beam elements and the lifting surface is divided into 14 elements in the chordwise direction and 160 elements in the spanwise direction to better capture the slipstream effect of the propeller. The calculation errors of the structural and aerodynamic models are within 1% compared with the situations where the element numbers are doubled.

According to the analysis framework of Section 2.1, the static aeroelastic vertical displacement of the wing at 4-degree angle-of-attack is shown in Fig. 13(a). Due to the acceleration effect of propeller slipstream, the maximum bending displacement of the wingtip increases by 5.8% compared with the clean wing. In the two cases with and without propellers, the wing bending displacement distribution calculated in this paper is consistent with that calculated by Werter [31] and Otsuka [32], and the static deformation error can be ignored.

Figure 13. Validation of the aeroelastic response: (a) static deformation, (b) dynamic response.

Using the static deformation above as the initial condition, a vertical force F = 100sin(ωt) N is applied at the wingtip to excite the vertical bending motion. ω is set to 2.24rad/s according to the first-order bending frequency of the structure. Figure 13(b) shows the vertical bending displacement response of the wingtip with and without propellers.

It can be seen that propeller slipstream increases the displacement amplitude. At the initial equilibrium position, the displacement difference caused by slipstream is about 0.25m. During the time history, the difference ΔMax of the maximum amplitude is about 0.17m and the difference ΔMin of the minimum amplitude is about 0.42m. Results indicate that the displacement difference caused by slipstream will gradually decrease/increase as the displacement response exceeds/falls below the equilibrium point. This is mainly because the stiffness of the geometric nonlinear structure will increase with deformation increasing. The displacement response trend calculated in this paper is consistent with that in Ref. [32], and the numerical errors of the two are mainly due to differences in the treatment of propeller aerodynamics.

5.0 Numerical results

In this section, the influence of distributed propellers on the dynamic response and aeroelastic stability of the clamped wing using the analysis framework above is investigated. The propeller-wing system studied is the same as that in Section 4.3. Different configurations are investigated including the clean wing, the wing with only thrust considered, and the wing with both propeller slipstream and thrust considered.

The configurations are calculated in the following three cases:

1. (1) For the response to bending excitation, a vertical step lumped force of 15N is applied on the wingtip to excite the bending dynamic motion.

2. (2) For the response to gust excitation, the wing is put under the gust with different frequencies. The effect of propeller speed is also investigated.

3. (3) For the stability study, the wing is put under the freestream with different velocities. The effect of propeller speed is also investigated.

Different initial conditions are applied for the three cases. For Case (1) and (2), as pointed out in Ref. [15], the dynamics of the highly flexible wing can be different for varying initial conditions because the structural stiffness will vary with the displacements. Therefore, the same initial displacements are artificially adopted for all configurations. For Case (3), zero deformation is taken as the initial displacement. For the used propellers with low speed and small mass, no gyroscopic effects are considered in the underlying structural model of all the cases investigated.

5.1 Response to bending excitation

Firstly, the influence of propellers on the response to bending excitation is explored. The freestream velocity is 25m/s, the air density is 0.0889kg/m³ and the angle-of-attack is 4°. A vertical step lumped force of 15N is applied on the wingtip elastic axis to excite the bending dynamic motion.

At the initial zero time, all configurations are under the same loads including aerodynamic force and propeller thrust of 1,900rpm. After that, thrust and slipstream are added or removed for different configurations during the time history. Figure 14 shows the time-domain response of the (a)vertical bending displacement and (b)torsion of the wingtip.

Figure 14. Response of the wingtip to bending excitation with a vertical step lumped force of 15N: (a) vertical bending displacement, (b) torsion.

As shown in Fig. 14(a), compared with the clean wing, propeller thrust reduces the maximum bending displacement by 11% and advances the phase by about 0.2 seconds. The distance between the propeller’s thrust axis and the wing’s elastic axis is zero. Therefore, this is mainly because the thrust is coupled with the nonlinear bending deformation of the highly flexible wing, resulting in the wing negative torsion and furtherly limiting the wing bending oscillation.

Compared with the configuration only considering thrust, propeller slipstream increases the maximum bending displacement by 3.2%. This is mainly because propeller slipstream increases the local velocity and thus increases the lift. Meanwhile, as shown in Fig. 14(b), propeller slipstream has little effect on the torsional response to bending excitation.

5.2 Response to gust excitation

The aerodynamic conditions such as the freestream velocity, the air density and the angle-of-attack are the same as those of Section 5.1. The initial state also remains unchanged. The three configurations are put under gust excitation.

The widely used 1-cos discrete gust model is selected to study the structural dynamic aeroelastic response of different configurations to gust excitation. The 1-cos discrete gust model can be expressed as follows, where w _g is the gust velocity amplitude and f _g is the gust frequency:

(34) \begin{equation}{w_{gust}} = \left\{ {\begin{array}{c@{\quad}c}{\dfrac{{{w_g}}}{2}\left( {1 - \cos 2\pi {f_g} \cdot t} \right)} {}& {0 \lt t \lt \dfrac{1}{{{f_g}}}}\\[5pt]0 & {}{t \gt \dfrac{1}{{{f_g}}}}\end{array}} \right.\end{equation}

In this section, efforts are devoted to cases with different propeller speed and different gust frequency. The specific parameter settings and their corresponding results are as follows:

1. (1) Set w _g = 1m/s, f _g = 1.5Hz. All propellers attached have the same speed and are set to 1,900, 2,000, 2,100, 2,200rpm in sequence. The vertical bending displacement response of the wingtip to gust excitation at different propeller speed is shown in Fig. 15.

   

   Figure 15. Vertical displacement response of the wingtip to gust excitation (w_g = 1m/s, f_g = 1.5Hz, and propeller speed = 1,900, 2,000, 2,100, 2,200rpm in sequence).

2. (2) Set all propeller speed to 2,000rpm and f _g to 0.5, 1.0, 1.5Hz in sequence. The vertical bending displacement response of the wingtip to gust excitation at different gust frequency is shown in Fig. 16.

   

   Figure 16. Vertical displacement response of the wingtip to gust excitation (propeller speed = 2,000rpm, w_g = 1m/s, and f_g = 0.5, 1.0, 1.5Hz in sequence).

The maximum vertical displacement response of the wingtip in different cases is shown in Fig. 17(a), and the response time corresponding to the maximum response is shown in Fig. 17(b).

Figure 17. (a) Maximum values and (b) the corresponding time of wingtip vertical displacement responses in different cases.

It can be seen that at the same gust frequency of 1.5Hz, with the increase of propeller speed, the negative torsion effect caused by propeller thrust is more evident. Compared with the clean wing, the maximum vertical displacement response of the wingtip decreases by 9.3%, 13.4%, 17.1% and 19.5%, respectively, and the corresponding time to the maximum response decreases by 5.3%, 10.6%, 25.5% and 41.3% respectively. Similarly, at the same propeller speed of 2,000 rpm, the maximum vertical displacement response and the corresponding response time decrease gradually with the increase of gust frequency. Propeller slipstream mainly increases the vertical displacement response and slightly increases the corresponding time to the maximum response (slightly delays the phase).

5.3 Influence on the aeroelastic stability

In order to study the influence of propeller thrust and slipstream on the aeroelastic stability of a highly flexible wing, the transient response at different freestream velocities is calculated. The equilibrium state at zero deformation is taken as the initial condition. The speed of all propellers is 2,000rpm and the angle-of-attack of the wing root is 2°. Figures 18 and 19 show the vertical bending and torsion response of the wingtip at the freestream velocity of 31.8 and 34.9m/s, respectively.

Figure 18. Transient response of the wingtip at the freestream velocity of 31.8m/s: (a) vertical displacement, (b) torsion.

Figure 19. Transient response of the wingtip at the freestream velocity of 34.9m/s: (a) vertical displacement, (b) torsion.

As shown in Fig. 18, at the freestream velocity of 31.8m/s, the bending and torsion responses of the clean wing oscillate with equal amplitude. However, under the same condition, the oscillation displacements of the wing converge due to the existence of propellers. As analysed in Section 5.1 above, propeller thrust limits the torsional deformation of the structure and the local angle-of-attack is reduced accordingly. Therefore, the energy obtained by the structure from aerodynamics is reduced, and the flutter speed is increased as a result.

As shown in Fig. 19, at the freestream velocity of 34.9m/s, the bending and torsion responses oscillate with equal amplitude when only considering the propeller thrust. Compared to the clean wing, propeller thrust increases the critical speed by 9.7%. When the influence of slipstream is introduced, the displacement response of the wing exhibits the characteristic of divergence compared to the response of the configuration with only thrust considered. This is mainly due to the fact that propeller slipstream increases the local flow velocity of the wing and leads to an increase in the aerodynamic energy obtained by the structure. Consequently, the flutter speed is decreased.

Let the calculation conditions remain unchanged. At different propeller speeds, the critical freestream velocities corresponding to the constant amplitude oscillation of the wingtip are shown in Fig. 20(a), and the oscillation frequencies are shown in Fig. 20(b). In the case of only considering thrust, as the propeller speed increases, the critical freestream velocity continuously increases while the oscillation frequency decreases. Compared with cases with only thrust considered, propeller slipstream reduces the critical velocity and increases the oscillation frequency. As previously commented, propeller slipstream brings an increase in local dynamic pressure, and causes a destabilising effect.

Figure 20. (a) Critical freestream velocity and (b) oscillation frequency at different propeller speeds.

6.0 Conclusions

In this paper a rapid nonlinear aeroelastic analysis framework had been established for the highly flexible wing considering the influence of distributed propellers attached along the spanwise direction. Nonlinear beam elements based on CR dynamics theory were applied for large deformation wing structure, and the local cylindrical coordinate systems of propellers in different position were derived to complete the coupling between the wing and the attached propellers. By taking advantage of the low-frequency dynamics of the high-aspect-ratio wing, the propeller wake was modeled as a quasi-steady skewed vortex cylinder without considering the wake updating process to improve calculation efficiency.

Validations for the isolated propeller and the propeller-wing unsteady system were conducted against published results and demonstrated the accuracy of the aerodynamic model as well as the rationality of the quasi-steady assumption of the propeller wake. The static and dynamic aeroelastic responses were also validated and results indicated that the framework proposed can effectively reflect the geometrically nonlinear characteristics of the highly flexible wing and the aerodynamic interaction characteristics of distributed propellers.

The framework was then employed to identify the effects of distributed propellers on the dynamic aeroelastic response and stability. From the cases, the following can be summarised:

1. (1) Large bending deformation caused thrust to produce wing negative torsion, which limited the aerodynamics and deformation. For the particular case, the maximum bending displacement was reduced by 11% and the phase was advanced by about 0.2 seconds compared with the clean wing. Additionally, the slipstream increased the aerodynamics and deformation by accelerating the local flow, thus increasing the maximum bending displacement by 3.2%. However, it had little effect on the torsional response to bending excitation.

2. (2) From the analysis of the response to gust excitation, it was observed that with the increase of propeller speed, the negative torsion effect of wing caused by thrust was more evident. Due to the presence of propellers, the maximum vertical displacement response and its corresponding response time were reduced by 9–42%, suggesting that distributed propellers can play a significant role in the gust alleviation of highly flexible aircraft. The propeller slipstream mainly increased the maximum values of the displacement response and slightly increased the corresponding time to the maximum response.

3. (3) The stability study showed that propeller thrust will limit the displacement divergence of the very flexible wing by reducing the structural torsional deformation and the local angle-of-attack. In addition, an improvement of flutter boundary was found with the increase of propeller speed. However, the propeller slipstream brought a destabilising effect as a result of the increment of dynamic pressure and local lift. Compared with the case with only thrust considered, the propeller slipstream decreased the critical speed by about 3%.

The thrust and slipstream of distributed propellers have a nonnegligible impact on the dynamic aerodynamic response of the highly flexible wing, making geometrically nonlinear effect more significant. In the initial design stage of such aircraft, these effects must be considered to avoid design deviations. The numerical investigations of the particular model point out the great potential of distributed propellers in gust alleviation and flutter suppression of highly flexible wings. In the subsequent research, the method of this paper can also be applied for the coupling design and optimisation between the very flexible wing and distributed propellers by adjusting parameters such as propeller speed, installation position, flexibility distribution, etc.