FLIGHT DYNAMIC-AEROELASTIC RESPONSE OF A HIGHLY FLEXIBLE AIRCRAFT WITH DISTRIBUTED PROPELLERS

International Forum on Aeroelasticity and Structural Dynamics IFASD 2024 17-21 June 2024, The Hague, The Netherlands FLIGHT DYNAMIC-AEROELASTIC RESPONSE OF A HIGHLY FLEXIBLE AIRCRAFT WITH DISTRIBUTED PROPELLERS Alberto Gallego Pozo1and Rauno Cavallaro1 1University Carlos III of Madrid Leganes, Madrid, Spain rauno.cav[email protected] Keywords: Aeroelasticity; highly-flexible Wings; Distributed Electric Propulsion; geometric nonlinearities; flight dynamic-aeroelastic coupling; whirl flutter. 1 INTRODUCTION 1.1 Aeroelastic Stability of Highly-Flexible Wings with Distributed Propellers: State of the Art Within the emergent electric aircraft market, Distributed Electric Propulsion (DEP) is a promising concept. It provides benefits in terms of aerodynamics and propulsive efficiency, noise reduction, and vehicle control [1]. Some of the new aircraft concepts applying DEP have high aspect ratio wings for improved efficiency. The light structure, in conjunction with the large propellers and motors mass and inertia, results in an enhanced structural flexibility of the wing, which can lead to an interaction between elastic deformations and rigid body flight dynamics that favors aeroelastic instabilities. This interaction may in turn be intensified by the gyroscopic effects induced by the propellers rotation, and the propellers aerodynamics. In effort [2], the aeroelasticity of a wing clamped at the root featuring distributed propellers is analyzed, and it is shown that the angular momentum of a large wing tip propellers has an important effect on the aeroelasticity of the wing. Contribution [3] explains that propeller whirl modes have an impact on the stability of the wing, after proposing a framework to study the aeroelastic stability of a wing clamped at the root featuring distributed propellers and linear beams to account for the structure. Effort [4] assessed the influence of whirl flutter on the design of the DEP Aircraft X-57 through multibody dynamic analyses on a semi-span clamped model. With the push towards more efficient aircraft, structural geometric nonlinearities are progressively becoming more relevant and need to be included in the design. These nonlinearities are generally a consequence of very large deflections but can also be driven early, at non-large deflections, by particular architectures, such as in Joined Wings cases [5]. Among the consequences of these nonlinearities at the aeroelastic level are the non-conservative prediction of static and dynamic instabilities and the relative change of mechanisms driving them [6]. The relevance of this phenomenon is exemplified by the Pazy wing test case [7], a benchmark model for geometrically nonlinear aeroelastic studies involving large deflections in low-speed flow, and by the vast body of literature about it. A thorough study of aeroelastic instabilities should consider possible coupling with flight dynamics responses. A phenomenon such as Body Freedom Flutter [8] is an example driven by 1 IFASD-2024-182 such interaction. On the other hand, even if no aeroelastic instability occurs, such coupling can impact flying qualities. The overall stability problem should be formulated considering a flying flexible body [9, 10]. In the current state of the art, there are studies on the aeroelastic assessment of distributed electric propulsion and on flight dynamic-aeroelastic coupling. However, no research has focused on highly flexible wings (and their relative nonlinearities) featuring distributed propulsion while accounting for the concurrent coupling with flight dynamics. 1.2 Challenges and Contributions The INDIGO (Integration and Digital Demonstration of Low-emission Aircraft Technologies and Airport Operations) project [11, 12], financed within the Horizon Europe programme, reunites academia, research centres and airports to identify the margins of improvement in airport Local Air Quality and Noise (LAQN) resulting from the introduction of a new non-conventional mid-range aircraft. The novel airframe features distributed propulsion based on hybrid electric/sustainable and conventional fuel powertrain and large aspect-ratio wing capable to fly quietly and in zero-to-low-emission mode (i.e. electric and SAF) at low altitudes near airports and resorts to conventional aviation fuel only when required, e.g., at higher altitudes or to recharge batteries during cruise. The concurrent integration of several propellers along a highly flexible wings increases the complexity of aeroelastic response for the following reasons: • Having typically lower natural frequencies, the aeroelastic modes have more room for interaction with the flight-dynamic ones. • With the wing being more flexible, the coupling with the classic whirl-flutter phenomenon is possibly tighter; the same holds for the coupling with the flight dynamics, where, at the propeller level, it is just a displacement and rotation of the mount point. • With the wing being more flexible, non-negligible differences arise when assuming a nondeformed reference condition or the real one in flight. This is true both at the aerodynamic level and the structural one (geometric nonlinearities). This study presents the derivation of the system of equations, incorporating all relevant coupled physics, and its implementation into a digital tool. The software’s capabilities are subsequently demonstrated using a synthetic baseline model of an aircraft with a highly flexible wing and DEP. Aeroelastic stability analyses are conducted using approaches of increasing complexity regarding the couplings and physical phenomena considered. The results are analyzed and discussed, emphasizing the differences in aeroelastic behavior when utilizing models that include specific couplings and physical effects. 2 THE COUPLED FLIGHT DYNAMIC-AEROELASTIC STABILITY MODEL This section presents the theoretical framework behind the coupled flight dynamic-aeroelastic stability module. Firstly, the perturbation equations of motion are derived for a general equilibrium condition. Then, the non-linear trim module to obtain the equilibrium condition is explained. Modal analysis is discussed to serve as appropriate basis for the analysis of the (small) displacements relative to the non-linear equilibrium. The aerodynamic models for the unsteady aerodynamics of both lifting surfaces and propellers are then described. Finally, the complete system of equations governing the flight dynamic-aeroelastic stability of the aircraft is derived and a solution procedure is suggested. 2 IFASD-2024-182 2.1 The Flight Dynamic-Aeroelastic Stability Equation The Stability Equation is derived starting from Lagrange’s equation: d dt(∂T ∂˙η)−∂T ∂η +∂U ∂η =∂(δW)) ∂(δη)(1) The formulation is oriented towards the integration of linear finite element models, which are the basis in aeroelasticity. The next assumptions are made for the derivation of the theoretical framework: •Assumption 1. The aircraft is discretized into finite elements using a lumped mass approach. Thus, the aircraft is discretized into points (k) with a mass (mk) and an inertia tensor (Jk). Each rotating mass (propeller) is then modelled as a point with a mass and an inertia tensor. •Assumption 2. Local translational and rotational elastic displacements with respect to an equilibrium condition are small; perturbation theory is used and linear elastic theory applies. •Assumption 3. Elastic displacements are described using an orthogonal mode shapes basis obtained from a free-free modal analysis of the aircraft. 2.1.1 Kinematics: position and velocity of a generic point of the aircraft Two reference frames are defined: an inertial frame (ΣI) attached to the Earth’s surface and a body frame (ΣB) attached to the aircraft. For a generic point of the aircraft: R=R0+r0+u(2) where: •R0≡position of the origin of ΣB •r0≡position of a differential of mass in the reference aircraft’s configuration •u≡variation of the position of the differential of mass with respect to the reference configuration due to deformations The next assumption is made: Assumption 4.dr0 dt B = 0 Let us differentiate between any generic point of the aircraft (denoted as with the sub-index k), all the points but the propeller points (sub-index i) and the points associated to each propeller (sub-index p). The linear velocities of points iand pare expressed in the same form, but their rotational velocities expressions differ: dRk dt I =v+ωB,I×(r0k+uk) + ˙ uk=v+ ˜ωB,I(r0k+uk) + ˙ uk Ωi=ωB,I+˙φi Ωp=ωB,I+˙φp+ωPH (3) where: 3 IFASD-2024-182 •˙ () = d dt()B •ωB,Iis the angular velocity of reference frame ΣBwith respect to ΣI. The adopted standard is that of the Tayt-Bryan Angles: ωB,IB =  p q r  =L  ˙ ϕ ˙ θ ˙ ψ  ;L=  1 0 −sinθ 0cosϕ cosθsinϕ 0−sinϕ cosθcosϕ  (4) •˙φiis the angular velocity of point idue to deformations; ˙φiis defined according to the Tayt-Bryan angles (each angle defines a rotational deformation angle at point i). These angles define a new reference frame ΣAiwith respect to ΣB. ˙φiAi =Li  ϕi θi ψi  ;Li=  1 0 −sinθi 0cosϕicosθisinϕi 0−sinϕicosθicosϕi  ;qroti=  ϕi θi ψi  (5) •˙φpis the angular velocity of point pdue to deformations and is defined, again, according to the Tayt-Bryan angles (each angle defines a rotational deformation angle at point p). These angles define a new reference frame ΣApwith respect to ΣB • The symbol ”˜” above a variable denotes the skew-symmetric matrix corresponding to a vector (for cross-products), such that: ˜ω=  0−ωzωy ωz0−ωx −ωyωx0  (6) •ωP,H is the angular velocity vector of the rotating mass at point p. 2.1.2 Kinetic energy The kinetic energy of the aircraft can be expressed as the sum between the translational kinetic energy and the rotational kinetic energy (T=EkinT+EkinR). EkinT=1 2X i ˙ RT i˙ Rimi+1 2X p ˙ RT p˙ Rpmp=1 2mvTv−X k mkvT(˜r0k+˜uk)ωB,I+X k mkvT˙uk− −1 2X k mkωB,IT(˜r0k+ ˜uk)(˜r0k+ ˜uk)ωB,I+X k mkωT B,I (˜r0k+ ˜uk)˙uk+1 2X k mk˙uT k˙uk(7) EkinR=1 2X i ΩiTJiΩi+1 2X p ΩpTJpΩp=1 2X k ωT B,IJkωB,I+X k ωB,ITJk˙φk+1 2X k ˙φT kJk˙φk+ +X p ωB,ITJpωPH +X p ˙φT pJpωPH +1 2X p ωPHTJpωPH (8) where Jkis the inertia tensor of point kand Gpis the gyroscopic matrix of the p−th propeller: Gp=ωPH   0Ixz −Ixy −Ixz 0Ix Ixy −Ix0  p (9) 4 IFASD-2024-182 2.1.3 Simplification of the kinetic energy expression Recalling assumptions 2and 3(small amplitude vibration around the equilibrium state), displacements can be written as a (truncated) superposition of (linearized about the reference system) of structural modes around the equilibrium state: u= [Φ]ηE(10) Mean axes [13] about the non-linear equilibrium state are used as body-reference frame, such that: • Origin of the body frame (mean axes frame) is located at the instantaneous center of mass: ZV rρdV = 0 −→ X k mk(r0+ud) = 0(11) • Internal linear momentum relative to elastic translations is zero: ZV ˙udρdV = 0 −→ X k mk˙udk=0(12) • Internal linearized angular momentum relative to elastic deformations is zero: ZV ˜r0˙udρdV = 0 −→ X k mk˜r0˙udk+X k Jk˙qrotdk=0(13) The previous two mean axes constraints are automatically satisfied since elastic displacements are approximated using the modal shapes of the free-free aircraft. ZV ˙udρdV = nE X i=1 dηi dt ZV ϕiρdV = 0 ZV ˜r0˙udρdV = nE X i=1 dηi dt ZV ˜r0ϕiρdV = 0 (14) Two further hypothesis which are common when analyzing the flight-dynamic-aeroelastic stability using mean axes are [14] [15] [16]: • Perturbation deformations and deformations rates are collinear: ZV ˜ud˙ud= 0 −→ X k mkωT B,I (˜r0k+ ˜uk)˙uk+X k ωB,IT[BRAk]JkLk˙qrotk= 0 (15) • The change in inertia due to perturbation deformations is negligible: 1 2X k ωB,IT(−mk(˜r0k+ ˜uk)(˜r0k+ ˜uk)+[BRAk]Jk[AkRB])ωB,I≈ ≈1 2X k ωB,IT(−mk˜r0k˜r0k+Jk)ωB,I=1 2ωB,ITJ0ωB,I(16) After these simplifications, the kinetic energy expression reads: T=1 2mvTv+1 2ωB,ITJ0ωB,I+1 2˙ηT EMEE ˙ηE+1 2X p ωPHTJpωPH+ +X p ωB,ITJpωPH −X p qrotp T[Gp]ωB,I+X p ˙qT rotpJpωPH +X p ˙qT rotp[Gp]qrotp (17) 5 IFASD-2024-182 2.1.4 Potential energy The potential energy of the aircraft is the sum of the gravitational potential energy Ugand the elastic strain energy Ue. Ug=−ZV ol (g·R)dm =−g·ROBm Uel =1 2ZV ol X ijkl CijklεijεkldV =1 2ηET[KEE]ηE (18) 2.1.5 Equations of motion The generalized coordinates that will be used are the position of the origin of the mean axes frame, the orientation angles of the mean axes frame relative the inertial one and the elastic modal coordinates describing displacements relative to the equilibrium. Applying Lagrange equations results in the following system of equations: m˙v +ωB,I0v+ ˜ωB,Iv0= ∂FB ∂p0 δp J0˙ωB,I+ ˜ωB,I0J0ωB,I+ ˜ωB,IJ0ωB,I0+X p [Gp]ωB,I+X p [Gp]˙qrotp+ +˜ωB,I0X p [Gp]qrotp= ∂MB ∂p0 δp MEE ¨ηE+KEEηE+X p [Φrotp]T[Gp][Φrotp]˙ηE+X p [Φrotp]T[Gp]ωB,I=∂δW ∂(δηE) (19) where vis now the perturbation velocity vector, v0the reference velocity vector, ωB,Ithe perturbation angular speed of the mean axis frame, ωB,I0the reference angular speed of the mean axis frame, and δpa vector containing all the perturbations that affect the forces and moments. For completeness, the kinematic equations need to be added to the previous system of equations. The inertial velocity of the aircraft is defined by the rate of displacement of its center of mass: V=vo+v=dvo dt I +dv dt I = ( ˙ XEiI+˙ YEjI+˙ ZEkI) + ( ˙xEiI+ ˙yEjI+ ˙zEkI) = = (UiB+VjB+WkB)+(uiB+vjB+wkB) (20) where the sub-index Idenotes inertial reference frame and Brefers to the mean axes frame. The orientation of the mean axes frame with respect to the inertial one is defined by the Euler angles (Φ = ϕ0+ϕ,Θ = Θ0+θ,Ψ = ψ0+ψ), according to the Tayt-Bryan Formalism. 6 IFASD-2024-182 2.2 Reference Condition: Trim 2.2.1 The aerodynamic problem The iterative nature of the non-linear flexible trim requires the use of a fast aerodynamic tool to evaluate the forces on the deflected configurations. Moreover, distributing rotary devices ahead of the wing results in complex aerodynamic interactions between the wake shed by the blades and the downstream surfaces. The selected compromise is to calculate the aerodynamics of the aircraft with DUST, a mid-fidelity tool developed to provide fast and reliable aerodynamic simulations of Vertical Take-Off and Landing (VTOL) aircraft configurations. The tool has been shown to provide reliable and fast predictions of the aerodynamic performance of unconventional VTOL aircraft [17–22]. The mathematical formulation behind this software relies on the Helmholtz’s decomposition of the velocity field, which allows to recast the aerodynamic problem as a combination of a boundary value problem for the potential part of the velocity and a mixed panels-vortex particles model for the free vorticity field in the flow. The lifting surfaces are modelled with surface panels, whereas propellers are modelled using lifting lines. Vorticity is shed from the trailing edge of both lifting lines and lifting surfaces and is then convected according to the local velocity, effectively considering the swirl imparted by the propellers on the flow impinging on the wing. The reader is referred to [23] for more details about the derivation. 2.2.2 The structural problem The structural in-house solver (pyBeam) [24] is based on a 6 DoF geometrically non-linear beam formulation. The Euler-Bernoulli beam kinematic assumption is considered. The equation governing the displacements of the structure in its discretized Finite Element (FE) form is: G(us) = fs−fint(us) = 0(21) where us,fsand fint are, respectively, the nodal generalized displacements, the external and internal load forces vector. The previous equation is solved by a Newton-Raphson method: Kus=−G(us)(22) where K=∂G(us) ∂usis the Jacobian/tangent matrix. 2.2.3 Splines and mesh deformation methods Aerodynamic and structural grids are generally non-coincident. A Moving Least Squares Algorithm is used to compute the spline matrix that relates structural and aerodynamic coordinates, displacements and forces [24]. Let xs∈RNsbe the coordinates of the structural nodes, and xa∈RNabe the coordinates of the aerodynamic moving grid, then it is possible to define a spline matrix S=S(us,ua), such that: ua=Sus fs=STfa (23) where uarepresents the displacements of the aerodynamic grid, usthe displacements of the structural grid, fsthe forces and moments on the structural nodes, and fathe forces and moments 7 IFASD-2024-182 on the aerodynamic nodes. The same procedure can be used to transfer displacements and forces among other relevant points of the grid, such as element centers. Within the present formulation, after solving for structural displacements, both the aerodynamic and structural grids can be updated by simply summing the displacements to the aerodynamic and structural grid coordinates of the previous iteration i−1: xai=xai−1+u∗ a−→ Da= 0 xsi=xsi−1+u∗ s−→ Ds= 0 (24) where a relaxation parameter αcan be applied to the boundary displacements to ensure stability of the method: u∗ a=αui a+ (1 −α)ui−1 a u∗ s=αui s+ (1 −α)ui−1 s (25) 2.2.4 The trim problem: Fluid-Structure Interaction (FSI) The equilibrium condition under consideration is steady-level flight. The equations that govern this equilibrium are: F(xa,c) = Lcosα −W Mycg =0;c=α δe(26) where αand δeare the Angle of Attack (AoA) and elevator (or any other control surface) deflection, respectively, and Mycg is the pitching moment with respect to the center of gravity. The equilibrium along the longitudinal direction (thrust equation) has been removed, since in a first approximation it is independent from the other two. Given the thrust that each propeller needs to produce, the necessary propeller’s collective pitch is first calculated using Blade Element Theory instead of the Vortex Particle Method to accelerate numerical computations, while retaining similar accuracy. Given an aerodynamic grid, the previous equation can be solved using a Good-Broyden’s method [25], which falls within the class of quasi-Newton methods. Due to the quasi linear relation between aerodynamic response, this method rapidly converges. The whole solution procedure to find the non-linear trim is described in the next image. 8 IFASD-2024-182 Figure 1: Non-linear trim workflow. 2.3 Modal Analysis After convergence of the trim procedure, the deflected shape is obtained. From the structural shape, it is possible to obtain the updated mass and stiffness matrices. A modal analysis is then performed to obtain the updated modes. Such shapes represent the modal basis to be used when formulating and resolving the perturbation aeroelastic equation. 2.4 Enhanced DLM for Potential Unsteady Aerodynamics The Doublet Lattice Method (DLM) is a vastly used method to compute the unsteady aerodynamics of lifting surfaces [26]. However, in its original formulation, it only calculates loads due to local pitching and plunging of the lifting surfaces, neglecting contributions due to in-plane motion and in-plane forces. This approximation provides good results on conventional wings. However, it may fail to provide accurate results on configurations where in-plane loads and the corresponding moments are important, like T-tails [27]. The in-plane forces can be relevant also when studying the free-free aircraft and important interactions between lateral-directional flight-dynamic and aeroelastic modes. In a high aspect-ratio wing, in-plane loads may cause non-negligible in-plane bending due to the higher flexibility. Furthermore, the gyroscopic motion of propellers tends to couple the yaw and pitch motion of the propeller’s axis (whirl mode). When the rotor’s mass is sufficiently large compared to the wing’s stiffness, this can lead to a coupling between in-plane bending, out-of-plane bending, and torsion. Hence, there may be a need to retain in-plane loads and motion in high aspect-ratio wings with distributed propellers. In the presence of large deflections the approximation of considering the aerodynamic forces about the undeformed configuration may lead to incorrect prediction. The effects of steady outof-plane bending and in-plane bending can be simply modelled by approximating the deformed lifting surface with several wing segments with different sweep and dihedral. However, the inclusion of angle of attack, twist and camber is more complex, due to the restriction that the velocity goes in the x−direction in the derivation of the pressure potential equation. Therefore a classic DLM cannot model these effects and modifications are required. To this aim, the general form of the non-penetration boundary condition (evaluated at the Control Point - CP - of each aerodynamic box j) needs to be derived [28]: Vj·nj= 0 −→ (U∞+u0j+u1jeiωt −iωhjeiωt)·(n0j+reiωt ×n0j)=0 (27) 9 IFASD-2024-182 where xR=u v w p q r xEyEzEϕ θ ψT To include the unsteady aerodynamics in the state-space system, lag states are defined: ik ik +βjV∞0 ˆceηR ηE=xlagj(6+nE)×1−→ ˙xlagj=L1L206×nE 0nE×12 InE×nExR ˙ηE−βjV∞0 ˆcxlagj(51) The new state vector is: x=    ˙xR ˙ηE ηE xlagj    (52) and the state-space system (for one lag state) is now: Ano lag 02(6+nE)×(6+nE) 0(6+nE)×2(6+nE)I6+nE    ˙xR ¨ηE ˙ηE ˙xlag1    = =               Bno lag       q∞0A2+1R 06×(6+nE) q∞0A2+1E 0nE×(6+nE)        L1L206×nE06×nE 0nE×12 InE×nE0nE  −β1V∞0 ˆcI6+nE                   xR ˙ηE ηE xlag1    (53) This system can be recasted into a classical eigenvalue-eigenvector analysis: (sI −A−1B)ξ=0(54) where the vector of state-space amplitudes ξcontains the amplitudes of rigid body states, aeroelastic states and artificial aerodynamic states. Modes tracking is required to identify the flight dynamic-aeroelastic eigenvalues and eigenvectors. The tracking mechanism is based on the iterative correlation of the eigenvectors for successive speeds. 3 APPLICATION TO A DEP AIRCRAFT WITH A HIGH ASPECT RATIO WINGS In this section, a synthetic testcase is proposed to demonstrate the capabilities of the presented framework. 3.1 Baseline aircraft The test model consists of a modified version of the NASA X-57, see Figure 7. The aerodynamic and structural models are based on a high aspect ratio wing with sweep angle, a rigid horizontal tailplane and a rigid vertical tailplane. The wing, HTP and VTP are rigidly connected to the center of mass of the aircraft. The structural FE model is a classic stick model, i.e., the structure is described by beams and rigid connections, as shown in Figure 8. The wing’s internal structure 16 IFASD-2024-182 Figure 7: NASA X57. is a single-cell aluminium wingbox. An added distributed mass is considered to account for non-structural mass. Figure 8: Structural grid. Figure 9: Propeller’s modelling in structural grid. The test case presents 2 large tip propellers and 4 smaller propellers distributed along the wing. The propellers are counter-rotating, with the propellers on the y−positive side rotating counterclockwise and those on the other side clockwise. The main properties of the test aircraft are Figure 10: Aerodynamic grid for trim procedure. collected in Tables 1, 2, 3. 17 IFASD-2024-182 Figure 11: DLM aerodynamic grid. Table 1: Wing, HTP and VTP geometry Wing HTP VTP span [m] 9.639 2.518 1.254 root chord [m] 0.756 0.783 1.773 tip chord [m] 0.334 0.783 0.699 Leading edge sweep [º] 9.879 0 48.469 Surface [m2]5.255 1.972 1.550 AR [-] 17.681 3.212 - Table 2: Aircraft mass properties Parameter Value aircraft mass [kg] 1174.82 wing mass (not including propellers) [kg] 144.70 total propellers mass [kg] 124.66 rest of aircraft mass [kg] 925.46 Xcg [m] 0.6821 Ycg [m] 0 Zcg [m] -0.5382 Ixx [kg·m2] 4386.923 Iyy [kg·m2] 6399.972 Izz [kg·m2] -115.952 18 IFASD-2024-182 Table 3: Propellers properties. Cruise propeller High-lift propeller Diameter [m] 1.024 0.387 root chord [m] 0.100 0.050 Number of blades 3 5 Airfoil MH117 MH114 Power [kW] 43 14.4 Rotational speed [rad/s] 2250 4548 Thrust [N] 578.98 220 (0 in cruise) Inclination angle with respect to wing [º] 0 0 rotating mass [kg] 8.12 2.00 non-rotating mass [kg] 36.77 5.44 Pitch stiffness [Nm] 12398.42 - Yaw stiffness [Nm] 14751.24 - 3.2 Test Cases The free-free stability analysis of the baseline aircraft is performed according to the cases listed in Table 4. The column Propeller indicates if gyroscopic and aerodynamic effects of Table 4: Test cases. Case Propellers Rigid-elastic coupling Reference condition DLM Reference aerodynamic shape Case 1 No No Undeformed Basic Undeformed Case 2 Yes No Undeformed Basic Undeformed Case 3 Yes Yes Undeformed Basic Undeformed Case 4 Yes No Deformed Basic Undeformed Case 5 Yes Yes Deformed Basic Undeformed Case 6 Yes Yes Deformed Enhanced Undeformed Case 7 Yes Yes Deformed Enhanced Deformed the propellers are included. The column Rigid-elastic coupling refers to including the flight dynamics/aeroelastic modes coupling. The column Reference condition specifies if the structural properties are evaluated on the undeformed configuration, or on the deflected trimmed one. The column DLM refers to the employment of the classic or enhanced versions of the DLM. And, the last column, Reference aerodynamic shape indicates whether the undeformed shape or the deflected one is used for the evaluation of the aerodynamic forces. Table 5: Cruise conditions V∞[m/s]h [m] ρ[kg/m3] 77.17 2438.4 0.9629 3.3 Results In this section, the results of the Flight Dynamic-Aeroelastic Stability analyses performed for all cases listed in Table 4 are presented, and a brief discussion is provided. 3.3.1 Normal modes Tables 6 and 7 summarize the results of the modal analysis on the aircraft when considering its undeflected (jig-shape) and deflected (in-flight) shapes. The propellers are not rotating, other19 IFASD-2024-182 wise, the modal analyses will results in out-of-phase mode shapes. Table 6: Natural modes of the unloaded aircraft (in its jig-shape configuration) Mode Frequency [Hz] 1S. First Bending 1.2973 1A. First Bending 2.1451 2S. First Torsion 5.238 2A. First Torsion 5.2453 3S. In-plane Bending 5.7043 4S. Second Bending 6.618 3A. Second Bending 6.792 4A. In-plane Bending 7.1448 5S. Tip propeller yaw 8.9689 5A. Tip propeller yaw 9.0702 Table 7: Natural modes of the aircraft in its in-flight deflected shape Mode Frequency [Hz] 1S. First Bending 1.2918 1A. First Bending 2.1814 2S. First Torsion 4.9998 2A. First Torsion 5.1424 3S. In-plane Bending 5.7137 3A. Second Bending 6.5152 4S. Second Bending 6.6025 4A. In-plane Bending 6.9997 5S. Tip propeller yaw 8.95 5A. Tip propeller yaw 9.0439 The letters (S) and (A) refer to symmetric and anti-symmetric modes; moreover, the shapes are described highlighting the most representative features: in fact, modes 2S/A to 4S/A are actually a combination of out-of-plane bending, in-plane bending, torsion and propeller pitch. A comparison between the two cases (undeflected and deflected structure) shows a general slight drop of the natural frequencies (with exception of mode 1A). A visual representation of the modes is provided in Section 5.1. It is interesting to observe that in-plane and outof plane motions are uncoupled in the planar case (jig-shape), but coupled in the deformed wing. In addition, the propeller pitch mode, which has a frequency of approximately 8 Hz when the mounting wing section is infinitely rigid, has a significant participation in the second bending and first torsion of both jig-shape and flight-shape cases. However, while it has a relevant participation in the in-plane bending of the jig-shape case, its participation in the inplane bending of the in-flight shape is negligible. 3.4 Stability analyses To support the discussion, several pictures and tables gathering the relevant data will be used. Figures 12 and 13 show the root loci of the first modes for Cases 1 and 2, respectively; together with Table 8 they will support the discussion on effects of wing-propeller interaction. Figures 14 is employed to highlight effects of flight dynamic-aeroelastic coupling in Cases 2 to 7. Tables 9 and 10 report the data of typical symmetric flutter occurrences. In particular, inspecting the root-loci for all Cases, it is inferred that flutter cases can be gathered in two groups, according to the instability frequency. These two groups are referred to with Type 1 (higher frequency) and Type 2 (lower frequency). Both Tables provide, for each flutter onset, the properties of the unstable aeroelastic mode in terms of participation of the modal basis. Hence, the real amplitude in the shape should be weighted according to the amplitude of modal base (which is mass normalized). 20 IFASD-2024-182 Figure 12: Case 1. Root Locus Figure 13: Case 2. Root Locus Table 8: Modes after introduction of gyroscopic effects at V∞= 0 m/s Modes 2S with Gyr. (| |,∠)Mode 2A with Gyr. (| |,∠)Mode 3S with Gyr. (| |,∠) 2S 0.94∠0◦2A 0.95∠0◦2S 0.56∠98◦ 3S 0.1∠86◦3A 0.04∠209◦3S 0.75∠0◦ 4S 0.05∠185◦4A 0.05∠92◦4S 0.06∠89◦ 5S 0.33∠91.04◦5A 0.32∠90◦5S 0.34∠177◦ Description S Backward whirl Description A Backward whirl Description S Backward whirl 21 IFASD-2024-182 3.4.1 Gyroscopic effects The frequencies on the pure imaginary axis, i.e., for no incoming flow, contain for Case 2 the effects of the inertial coupling with the rotating propellers. Table 8 shows that the modes 2S, 2A, and 3S for Case 2, i.e., including the gyroscopic effects, are pretty different than the ones of Case 1. For example, the new 2S mode, is now a superposition of the original modes 2S, 3S, 4S and 5S, each with a different phase; this mode now features a strong backward whirl component. Overall, the inclusion of gyroscopic effects couple in-plane bending, propeller pitch, yaw, and torsion, resulting in effective whirl modes of the propeller. As a consequence, the frequencies of modes 2S and 3S largely decrease. 3.4.2 Wing-propellers aeroelastic coupling Inspection of the root loci for Case 1 and 2 (Figures 12 and 13) depicts a complex picture. In Case 1, three instability occurrences are detected for the first 4 modes. For the symmetric case, mode 2S and 3S become unstable, even though mode 2S shows a hump-mode flutter. The frequencies of the flutter are pretty close, around 5.7 Hz, suggesting possible complex interactions. For the anti-symmetric case, mode 2A is the one becoming unstable. Flutter frequency is also similar to the one of the symmetric case. In Case 2, i.e., introducing the wing-propellers coupling, the flutter picture is pretty different. Both modes 2S and 3S, featuring backward-whirls, become unstable, but at pretty different flutter frequencies ( 4.2 and 5.9 Hz, approximately). For the 2S instability, the flutter speed is reduced of 10%. The flutter mechanism changes from being a purely torsion-bending coupling to being a combination of torsion-bending and an effective backward whirl of the tip propeller, including the propeller yaw (see Table 10). For the instability of mode 3S, a similar flutter mechanism is observed. However, the flutter speed is largely increased and in-plane motion becomes more relevant, (see Table 9). Mode 2A, featuring a backward-whirl, also flutters slightly over 4 Hz. 3.4.3 Flight dynamic and aeroelastic coupling Attention is now drawn to the differences among Cases 2 to 7 concerning the flight dynamic modes, i.e., the eigenvectors that resemble the classic flight-dynamic modes of a rigid aircraft. Figure 14 shows the damping and frequency of the short period and the dutch roll modes of the flexible aircraft. 22 IFASD-2024-182 Figure 14: Damping and frequency response of flight dynamic modes of the flexible aircraft. The frequency evolution of these modes is observed to be nearly independent of the fidelity of the approach. For example, Case 2, which discards the effects of flexibility on flight dynamic response, already predicts the frequency reasonably well. This is not the case for the damping of the short-period mode. Starting from Case 3, in which the flight dynamic/aeroelastic coupling is taken into account, a decrease in damping is observed. This is possibly a consequence of the typical coupling between the short-period mode and the first symmetric bending mode. Additionally, the contributions introduced by the EDLM also reduce the damping of the shortperiod mode (Case 5 vs. Cases 6 7). With regards to the damping of the Dutch roll, the reader may notice a decrease at low speeds for Case 7, leading to instability. For the other cases, damping is only marginally changing. A possible explanation is that a deflected wing is equivalent to a wing with more dihedral, resulting in larger roll damping. This, in turn, increases the difference in magnitude between roll and yaw damping, thereby reducing the overall stability of the Dutch roll. 3.4.4 Flutter Type 1 From this point on, only the symmetric modes are considered for the sake of clarity and conciseness. As mentioned above, from Case 2 to Case 7, all flutter occurrences have been noted to happen around two frequency ranges. Type 1 includes all flutter occurrences with a frequency of about 5.5-5.9 Hz. The mode that loses stability is Mode 3S, which, at zero wind speed, primarily features in-plane bending, torsion, and tip-propeller yaw, resembling a backward-whirl mode. Figure 15 shows mode 3S damping and frequency vs speed. From inspection, it can be inferred that updating the structural part of the aeroelastic system with the real stiffness and inertial distribution (geometric nonlinearities), i.e., from Case 4 onward, has a significant impact on the flutter onset. The frequency of the aeroelastic mode increases by approximately 0.2 Hz or more, while damping is reduced, leading to an earlier onset of flutter (flutter speed drops by about 100 m/s). From Case 4, the flutter mechanism also appears to change, with the unstable aeroelastic mode showing an increased participation of Mode 2S 23 IFASD-2024-182 (first torsion) and an almost negligible participation of Mode 4S (second bending). Moreover, looking at the phases of the 2S, 3S and 5S modes suggests that flutter mode is governed by an effective backward whirl mode of the propeller’s hub (note the differences in the phases of 2S and 5S between Cases 2 and 4). It should be noted, however, that the natural modes evaluated on the deflected shape do not closely resemble those evaluated on the jig shape, see Section 5.3. The reader is referred to Section 5.5 for the modes of deformed shape including gyroscopic effects. Further analyses are ongoing to gain a deeper understanding of the results. Figure 15: Damping and frequency of aeroelastic mode 3S. Table 9: Flutter type 1. Symmetric case. Flutter mechanism, frequency, and speed. Flutter Type 1 Natural Mode Case 1 Case 2 Case 3 Case 4 Case 5 Case 6 Case 7 w0 0 1.38∠15◦0 0.05∠265◦0.07∠280◦0.05∠282◦ q0 0 0.18∠158◦0 0.04∠271◦0.04∠265◦0.03∠259◦ 1S. 1st Bending (S) 0.10∠181◦0.25∠255◦0.26∠255◦0 0 0 0 2S. 1st Torsion (S) 1.00∠0◦0.53∠124◦0.51∠127◦0.86∠35◦0.88∠35◦0.86∠35◦0.86∠34◦ 3S. In-Plane Bending (S) 0.30∠8◦1.00∠0◦1.00∠0◦1.00∠0◦1.00∠0◦1.00∠0◦1.00∠0◦ 4S. 2nd Bending (S) 0.43∠168◦0.93∠234◦1.00∠234◦0.06∠206◦0.06∠206◦0.03∠206◦0.03∠200◦ 5S. Propeller yaw (S) 0 0.60∠169◦0.62∠169◦0.37∠221◦0.38∠221◦0.36∠221◦0.35∠220◦ Flutter Mechanism 2nd Bending - 1st torsion - In-plane bending (S) In-plane bending - 1st Torsion - Prop.Yaw (S) (Equivalent to a backward whirl of tip propeller) Flutter Frequency [Hz] 5.75 5.89 5.89 5.65 5.66 5.65 5.65 Flutter Speed [m/s] 92.75 152.25 158 51.74 50.75 59.25 62.25 3.4.5 Flutter Type 2 Type 2 includes all flutter occurrences with a frequency of about 4.2-4.3 Hz. The mode that loses stability is the 2S, which, at zero wind speed, primarily features torsion and tip-propeller yaw, resembling a back-whirl mode. Figure 16 shows mode 2S damping and frequency vs speed. 24 IFASD-2024-182 Focusing on the trend of damping, three main jumps are noticed. The first one, from Case 3 to 4, is a stabilizing effect induced by updating the structural part of the aeroelastic system with the real stiffness and inertial distribution (geometric nonlinearities). The second jump, again stabilizing is observed when including the flight dynamic-aeroelastic coupling (from Case 4 to 5, but also from Case 2 to 3). The third jump, promoting instability, is the introduction of the aerodynamic effects typically neglected in the classic DLM (from Case 5 to Cases 6 and 7). The flutter mechanism changes between Case 2 & 3 and Cases 6 & 7. Now, mode 1S (first bending) and 4S (second bending) participate in the unstable aeroelastic mode. Figure 16: Damping and frequency of aeroelastic mode 2S. Table 10: Flutter type 2, symmetric case. Flutter mechanism, frequency, and speed. Flutter Type 2 Natural Mode Case 1 Case 2 Case 3 Case 4 Case 5 Case 6 Case 7 w0 0 0.69∠296◦- - 2.64∠275◦2.12∠272◦ q0 0 0.08∠104◦- - 0.20∠135◦0.21∠133◦ 1S. 1st Bending (S) 0.09∠181◦0.17∠185◦0.19∠186◦- - 0.92∠173◦0.79∠172◦ 2S. 1st Torsion (S) 1.00∠0◦1.00∠0◦1.00∠0◦- - 1.00∠0◦1.00∠0◦ 3S. In-Plane Bending (S) 0.58∠159◦0.13∠91◦0.13∠92◦- - 0.33∠94◦0.31∠89◦ 4S. 2nd Bending (S) 0.36∠169◦0.19∠184◦0.19∠185◦- - 0.85∠352◦0.71∠351◦ 5S. Propeller yaw (S) 0 0.40∠93◦0.40∠93◦- - 0.69∠81◦0.65∠82◦ Flutter Mechanism 2nd Bending - 1st Torsion - In-plane Bending (S) 1st Torsion - Prop. Yaw (S) (Equivalent to a backward whirl of tip propeller) - - Out-of-plane Bending - 1st Torsion Flutter Frequency [Hz] 5.67 4.24 4.24 - - 4.28 4.26 Flutter Speed [m/s] 88.75 81.75 86.25 - - 151.25 143.75 4 CONCLUSIONS In this paper, the formulation of a framework for the aeroelastic stability assessment of highlyflexible configurations featuring distributed electric propulsion is presented. The approach considers the aircraft free in the air, thus retaining the flight dynamic-aeroelastic coupling. Additionally, it takes into account the aeroelastic effects induced by the propeller as well as the effects of large deflections. To this aim, an ad-hoc enhanced version of the DLM has been developed and integrated. Moreover, a trim procedure has been established to find the reference 25 IFASD-2024-182 Figure 28: Mode 2S. First Symmetric Torsion Figure 29: Mode 3S. First Symmetric In-Plane Bending 32 IFASD-2024-182 Figure 30: Mode 4S. Second Symmetric Bending Figure 31: Mode 5S. Symmetric Propeller yaw 33 IFASD-2024-182 5.4 Antisymmetric modes Figure 32: Mode 1A. First Antisymmetric Bending Figure 33: Mode 2A. First Antisymmetric Torsion 34 IFASD-2024-182 Figure 34: Mode 3A. Second Antisymmetric Bending Figure 35: Mode 4A. First Antisymmetric In-Plane Bending 35 IFASD-2024-182 Figure 36: Mode 5A. Antisymmetric Propeller yaw 5.5 Gyroscopic modes for the in-flight aircraft Table 11: Modes of the deformed aircraft after introduction of gyroscopic effects at V∞= 0 m/s Natural mode Gyr mode 1 Natural Mode Gyr mode 2 Natural Mode Gyr mode 3 2S 0.88∠0◦2A 0.94∠0◦2S 0.64∠35◦ 3S 0.28∠129◦3A 0.06∠272◦3S 0.72∠0◦ 4S 0.08∠345◦4A 0.05∠318◦4S 0.08∠213◦ 5S 0.38∠78◦5A 0.33∠87◦5S 0.27∠223◦ Description S Backward whirl Description A Backward whirl Description S Backward whirl Frequency [Hz] 4.14 Frequency [Hz] 4.19 Frequency [Hz] 5.65 6 ACKNOWLEDGEMENTS The activities described in this paper have been carried out under the project INDIGO (Integration and Digital Demonstration of Low-emission Aircraft Technologies and Airport Operations), coordinated by R. Cavallaro from Universidad Carlos III de Madrid. INDIGO project [11, 12] has received funding from the European Climate, Infrastructure and Environment Executive Agency (CINEA) under the Horizon Europe programme under grant agreement No 101096055. 7 REFERENCES [1] Kim, H. D., Perry, A. T., Ansell, et al. A review of distributed electric propulsion concepts for air vehicle technology. [2] Amoozgar, M., Friswell, M. I., Fazelzadeh, S. A., et al. (2021). Aeroelastic stability analysis of electric aircraft wings with distributed electric propulsors. Aerospace, 8. ISSN 22264310. doi:10.3390/aerospace8040100. [3] B¨ ohnisch, N., Braun, C., Koschel, S., et al. (2022). Whirl flutter for distributed propulsion systems on a flexible wing. American Institute of Aeronautics and Astronautics Inc, AIAA. ISBN 9781624106316. doi:10.2514/6.2022-1755. 36 IFASD-2024-182 [4] Hoover, C. B., Shen, J., Kreshock, A. R., et al. Whirl flutter stability and its influence on the design of the distributed electric propeller aircraft x-57. [5] Cavallaro, R. and Demasi, L. (2016). Challenges, ideas, and innovations of joined-wing configurations: A concept from the past, an opportunity for the future. Progress in Aerospace Sciences, 87, 1–93. ISSN 0376-0421. doi:https://doi.org/10.1016/j.paerosci. 2016.07.002. [6] Gray, A. C., Riso, C., Jonsson, E., et al. (2023). High-fidelity aerostructural optimization with a geometrically nonlinear flutter constraint. AIAA Journal, 61(6), 2430–2443. doi: 10.2514/1.J062127. [7] Avin, O., Raveh, D. E., Drachinsky, A., et al. (2022). Experimental aeroelastic benchmark of a very flexible wing. AIAA Journal, 60(3), 1745–1768. doi:10.2514/1.J060621. [8] Love, M. H., Zink, P. S., Wieselmann, P. A., et al. (2005). Body freedom flutter of high aspect ratio flying wings. AIAA 2005-1947. 46th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Material Conference, Austin, Texas. doi:10.2514/6. 2005-1947. [9] Milne, R. D. (1962). Dynamics of the Deformable Airplane. Her Majesty’s Stationary Office, (3345). [10] Bombardieri, R., Cavallaro, R., Castellanos, R., et al. (2021). On the dynamic fluid–structure stability response of an innovative airplane configuration. Journal of Fluids and Structures, 105, 103347. ISSN 0889-9746. doi:https://doi.org/10.1016/j.jfluidstructs. 2021.103347. [11] Innovative aircraft and propulsion technologies to improve air quality near airports. doi: https://doi.org/10.3030/101096055. [12] Innovative aircraft and propulsion technologies to improve air quality near airports. https://indigo-sustainableaviation.eu/. [13] Canavin, J. R. and Likins, P. W. (1977). Floating reference frames for flexible spacecraft. Journal of Spacecraft and Rockets, 14, 724–732. ISSN 00224650. doi:10.2514/3.57256. [14] Waszak, M. R. and Schmidt, D. K. (1988). Flight dynamics of aeroelastic vehicles. Journal of Aircraft, 25, 563–571. ISSN 00218669. doi:10.2514/3.45623. [15] Schmidt, D. K. D. K. (2012). Modern flight dynamics. McGraw-Hill. ISBN 007339811X. [16] Carey S. Buttrill, T. A. Z. and Arbuckle, P. D. Nonlinear simulation of a flexible aircraft i n maneuvering flight. [17] Tugnoli, M., Montagnani, D., Syal, M., et al. (2021). Mid-fidelity approach to aerodynamic simulations of unconventional vtol aircraft configurations. Aerospace Science and Technology, 115. ISSN 12709638. doi:10.1016/j.ast.2021.106804. [18] Palazzi, M. (2020). Mid-fidelity approach to tiltrotor aerodynamics. [19] Montagnani, D., Tugnoli, M., Fonte, F., et al. (2019). Mid-fidelity analysis of unsteady interactional aerodynamics of complex vtol configurations. 37 IFASD-2024-182 [20] Savino, A., Cocco, A., Zanotti, A., et al. (2021). Coupling mid-fidelity aerodynamics and multibody dynamics for the aeroelastic analysis of rotary-wing vehicles. Energies, 14. ISSN 19961073. doi:10.3390/en14216979. [21] Perdolt, D., Thiele, M., Milz, D., et al. (2021). Comparison of multi-fidelity rotor analysis tools for transitional and low speed flight regimes. doi:10.25967/550128. [22] Fern´ andez, S. C. (2023). Mid-fidelity approach to propeller-wing interactional aerodynamics. [23] Uhlman, J. S. An integral equation formulation of the equations of motion of an incompressible fluid. [24] Bombardieri, R., Cavallaro, R., Sanchez, R., et al. (2021). Wing shape optimization assisted by algorithmic differentiation. Structural and Multidisciplinary Optimization. [25] Broyden, C. G. (2000). On the discovery of the “good broyden” method. Springer-Verlag. Math. Program., Ser. B, 87, 209–213. doi:10.1007/s101070000155. [26] Albano, Y. E., Northrop, W. P. R., Hawthorne, N., et al. No. 68-73-a doublet lattice method for calculating lift distributions on oscillating surfaces in subsonic flows meeting a doublet lattice method for calculating lift distributions on oscillating surfaces in subsonic flows. [27] Zy, L. H. V. and Mathews, E. H. (2011). Aeroelastic analysis of t-tails using an enhanced doublet lattice method. Journal of Aircraft, 48, 823–831. ISSN 15333868. doi:10.2514/ 1.C001000. [28] Zyl, L. H. V. (2008). Unsteady panel method for complex configurations including wake modeling. Journal of Aircraft, 45, 276–285. ISSN 15333868. doi:10.2514/1.29267. [29] Kier, T. M. and Looye, G. H. N. Unifying manoeuvre and gust loads analysis models. [30] Bombardieri, R. (2021). Aerostructural optimization and aeroelasticity of new generation aircraft. Ph.D. thesis. [31] Roger, K. L. (1977). Airplane math modeling methods for active control design. AGARDCP-228, 4–1. [32] Baldelli, D. H., Chen, P. C., and Panza, J. (2006). Unified aeroelastic and flight dynamic formulation via rational function approximations. Journal of Aircraft, Vol. 43(No. 3), pp. 763–772. doi:http://dx.doi.org/10.2514/1.16620. [33] Dettman, J. (1969). Mathematical Methods in Physics and Engineering. Dover Books on Engineering. Dover. ISBN 9780486656496. [34] Reed, W. H. (1967). Review of propeller-rotor whirl flutter. [35] Rose, T. and Rodden, B. (1989). Propeller/nacelle whirl flutter addition on msc nastran. [36] III, W. and Bland, S. (1961). An analytical treatment of aircraft propeller precession instability. 38 IFASD-2024-182 COPYRIGHT STATEMENT The authors confirm that they, and/or their company or organisation, hold copyright on all of the original material included in this paper. The authors also confirm that they have obtained permission from the copyright holder of any third-party material included in this paper to publish it as part of their paper. The authors confirm that they give permission, or have obtained permission from the copyright holder of this paper, for the publication and public distribution of this paper as part of the IFASD 2024 proceedings or as individual off-prints from the proceedings. 39