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 [Reference Otsuka, Del and Palacios24]. 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 [Reference Branlard and Gaunaa25]. 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 [Reference Patil and Hodges13–Reference Hodges, Patil and Chae16]. 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 [Reference Crisfield, Galvanetto and Jelenic26].

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.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 [Reference Veldhuis27] 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 [Reference Branlard and Gaunaa25] 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 (Reference Crisfield, Galvanetto and Jelenic26) 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. [Reference Crisfield, Galvanetto and Jelenic26].

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 (Reference Teixeira and Cesnik22) 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. [Reference Xie and Yang28].

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.1 Validation of the isolated propeller

As a first check for the implementation of propeller aerodynamics, the comparison of thrust with results from Ref. [Reference Fan, Zhou, Zhu, Wang and Wang29] 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. [Reference Fan, Zhou, Zhu, Wang and Wang29], 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. [Reference Veldhuis27], 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. [Reference Xue30], 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. [Reference Otsuka, Del and Palacios24]. 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. [Reference Otsuka, Del and Palacios24], 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 [Reference Werter and Breuker31, Reference Otsuka, Carre and Palacios32] 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 [Reference Werter and Breuker31] and Otsuka [Reference Otsuka, Carre and Palacios32], 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. [Reference Otsuka, Carre and Palacios32], and the numerical errors of the two are mainly due to differences in the treatment of propeller aerodynamics.

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.