T-tail flutter simulations with regard to quadratic mode shape components

It is known that the dynamic aeroelastic stability of T-tails is dependent on the steady aerodynamic forces at aircraft trim condition. Accounting for this dependency in the flutter solution process involves correction methods for doublet lattice method (DLM) unsteady aerodynamics, enhanced DLM algorithms, unsteady vortex lattice methods (UVLM), or the use of CFD. However, the aerodynamic improvements along with a commonly applied modal approach with linear displacements results in spurious stiffness terms, which distort the flutter velocity prediction. Hence, a higher order structural approach with quadratic mode shape components is required for accurate flutter velocity prediction of T-tails. For the study of the effects of quadratic mode shape components on T-tail flutter, a generic tail configuration without sweep and taper is used. Euler based CFD simulations are applied involving a linearized frequency domain (LFD) approach to determine the generalized aerodynamic forces. These forces are obtained based on steady CFD computations at varying horizontal tail plane (HTP) incidence angles. The quadratic mode shape components of the fundamental structural modes for the vertical tail plane (VTP), i.e., out-of-plane bending and torsion, are received from nonlinear as well as linear finite element analyses. Modal coupling resulting solely from the extended modal representation of the structure and its influence on T-tail flutter is studied. The g-method is applied to solve for the flutter velocities and corresponding flutter mode shapes. The impact of the quadratic mode shape components is visualized in terms of flutter velocities in dependency of the HTP incidence angle and the static aerodynamic HTP loading.


Introduction
Environmental as well as economical aspects motivate aircraft jet engine manufacturers to continuously improve their jet engine's efficiency. An increase in bypass ratio is most often a feasible way to reduce fuel consumption and noise emission. The associated larger engine diameters, however, pose a challenge for the engine integration under an aircraft's wing due to the limited ground clearance. A solution to this, most commonly found in business aviation, is the placement of the engines at the rear fuselage in combination with a T-tail empennage.
For dynamic aeroelastic considerations, T-tails are in a certain way unusual compared to conventional tails in that the stability strongly depends on the static aerodynamic loading and unsteady aerodynamic forces induced by inplane motion, which can usually be neglected in flutter simulations. A typical T-tail flutter mechanism can be described as a coupling between elastic VTP and HTP deformations, usually involving out-of-plane bending and torsion of the VTP. As the HTP performs roll motion when the VTP undergoes out-of-plane bending, the steady lift force vector is tilted in lateral direction, inducing a lateral force component. For VTP torsion, the HTP yaws, which results in asymmetrical aerodynamic span loading and a corresponding rolling moment. An assessment of T-tail flutter stability therefore requires correction of DLM unsteady aerodynamics [1][2][3], the use of enhanced DLM algorithms [4], UVLM [5], or CFD approaches [6][7][8][9][10]. Van Zyl pointed out, however, that the consideration of the additional aerodynamic terms alone 1 3 can actually reduce the flutter velocity prediction accuracy when a modal approach is applied with a linear displacement model. He suggests the use of an extended modal approach with quadratic mode shape components for a kinematically reasonable deformation model [4,[11][12][13].
The concept of quadratic mode shape components has been introduced by Segalman and Dohrmann for rotating flexible structures in [14][15][16][17][18]. In their works, they propose a series of nonlinear static finite element analyses to obtain the quadratic mode shape components. Ritter uses this idea for simulating the nonlinear aeroelastic gust responses of free-flying aircraft in [19]. Van Zyl seized upon the subject of quadratic mode shape components for T-tail flutter studies and proposes a method for computing the quadratic mode shape components based on linear finite element analyses [12]. This approach has been applied most recently to T-tail simulations by Murua in combination with potential flow theory based aerodynamics [5]. The approach based on linear finite element analyses is also applied by Farao et al. for the simulation of nonlinear gust response of a full aircraft model [20].
Previous work of the author regarding T-tail flutter addressed an assessment of the gain in flutter velocity prediction accuracy obtained by augmenting DLM unsteady aerodynamic forces by means of a strip theory approach in comparison to aerodynamic forces from LFD simulations [21]. While quadratic mode shape components have been neglected so far, this paper focuses on T-tail flutter simulations with quadratic mode shape components and CFD based unsteady aerodynamic forces. The CFD approach, while being computationally more expensive compared to a DLM approach, inherently captures the aerodynamic forces to the full extend necessary for T-tail flutter assessment and, hence, does not require correction for lateral aerodynamic forces induced by roll motion of the HTP as well as unsteady aerodynamic forced induced by inplane motion of the HTP. The quadratic mode shape components are obtained by a series of nonlinear as well as linear static finite element analyses. In contrast to the cited works regarding T-tail flutter, both, uncoupled as well as coupled quadratic mode shape components, are considered for the flutter assessment. Here, the focus is put on the coupling due to the quadratic mode shape components resulting from the first VTP out-of-plane bending and first VTP torsion.

Modal analysis
The modal analysis for obtaining the structural eigenvalues and eigenvectors is performed using the solution sequence SOL 103 provided by MSC.NASTRAN [22]. By applying the Lanczos method, the real eigenvalue problem (Eq. (1)) is solved for the structural eigenmodes and eigenfrequencies. K denotes the structural stiffness matrix and M the structural mass matrix in the physical reference frame. i is the angular eigenfrequency and i the corresponding real eigenvector.
Combining the eigenvectors i into a mode shape matrix with and m being the number of mode shapes, leads to the generalized mass and stiffness matrices in the modal reference frame, viz.

g-Method flutter equation
The generalized flutter equation in the Laplace domain can be expressed in terms of the generalized mass matrix M , the generalized stiffness matrix K , the generalized aerodynamic forces matrix Q(p, M ∞ ) , and the vector of generalized coordinates q as [23] where p = k + jk = g + jk is the non-dimensional Laplace parameter, j the imaginary unit, the decay rate coefficient, q ∞ the dynamic pressure, and the reduced frequency parameter. c∕2 and V denote the reference length and velocity, respectively, and M ∞ the Mach number.
The generalized aerodynamic forces are usually valid only for simple harmonic motion (i.e., p = jk ) and, thus, only at flutter onset. Below and above this point, the g-method [24] assumes a formulation of the generalized aerodynamic forces for small damping values g according to where Q � (jk, M ∞ ) = dQ(jk, M ∞ )∕d(jk) . This assumption leads to the g-method flutter equation, viz.
A solution of Eq. (8) is obtained by rewriting it in terms of an eigenvalue problem with x being the eigenvector and As only real valued damping factors are of interest, a reduced frequency sweep searches for a sign change in the imaginary part of g with increasing reduced frequency. At the sign change, the flutter frequency and damping values are computed based on a linearly interpolated reduced frequency.

Linearized frequency domain solver
The LFD solver implemented in the DLR TAU code [25][26][27][28][29][30] is used for small perturbation simulations w.r.t. a (nonlinear) reference state with linearized Euler and Navier-Stokes equations. The finite volume discretized RANS equations are written in form of an integral over the control volume surface and an integral over the control volume with W = (1, u, v, w, E, ) T as the vector of fluid unknowns, x as a vector of grid coordinates, and R as the residual function. f c ⋅ n is the convective and f v ⋅ n the viscous flux, while n denotes the surface normal vector and Q the turbulent source vector. In semi-discrete form and with as a diagonal matrix containing the cell volumes, this equation reads For the linearization, the amplitudes of the grid motion x and the vector of fluid unknowns W are assumed to be small and, hence, can be expressed as a Taylor series expansion truncated after the first-order term, namely By assuming a converged steady reference state with Eq. (15) becomes Transformation of Eq. (19) in the frequency domain is achieved by expressing the perturbation vectors of fluid unknowns and grid displacements as Fourier series truncated after the first harmonic, viz.
These considerations yield the linear system of equations The residual Jacobian R∕ W on the left hand side of Eq. (23) is computed analytically, while the right hand side is solved for by applying a central difference scheme with a deformed mesh at positive as well as negative grid deformation amplitude. For the propagation of the surface mesh deformation into the grid volume, Radial basis functions are used [31,32].
In comparison with full nonlinear CFD computations, the increase in computational performance comes with a loss of accuracy when nonlinear aerodynamic effects are induced by the small perturbations.

Extended modal approach
The physical displacements of a structural grid point are formulated in terms of the linear mode shapes i and the corresponding quadratic mode shape components g ij , viz.
where q i,j are the ith and the jth modal degree of freedom, respectively. The first term in Eq. (24) may be obtained by conventional modal analysis (see section 2.1). For the second term, two methods will be employed in the present work, where a distinction is being made between a nonlinear and a linear static finite element approach.

Obtaining quadratic mode shape components from nonlinear static analysis
To evaluate the quadratic mode shape components, a series of modal force fields is defined. The factors s i,j a,b are used to scale the linear modes i and j and are centered around zero. For small deformations and i = j , the deformations are equivalent to the (scaled) mode shape i. In analogy to [19], the number of scaling factors will be referred to as the stencil size . Applying these force fields in terms of a series of nonlinear static analyses results in a vector of nonlinear solutions N with 2 rows. A quartic polynomial is then fitted through the resulting nonlinear deformations, which leads to a linear system of 2 equations and 14 unknowns, viz.
with A solution of the usually over-determined linear system of Eq. (26) may be found in terms of a least squares approach. The solution vector u ij then contains 14 rows from which the quadratic mode shape components can be extracted, as indicated by the columns with quadratic components in the P ij -matrix where x ij is the vector of physical displacements.

Obtaining quadratic mode shape components from linear static analysis
In a linear normal modes analysis, the rotation of an element is expressed in terms of the linear displacements of its nodes, which may lead to element stretching. For a finite structural 1-D element k, hereinafter referred to as line element k, in linear mode i, this stretching can be expressed as [12] where l k is the line element vector and is the element rotation vector. which is solved for the uncoupled quadratic mode shape component The Lagrange multiplier is used for imposing an orthogonality condition on the quadratic mode shape component w.r.t. the stiffness matrix For the coupled quadratic mode shape components g ij , a similar procedure is followed with an artificial mode shape consisting of a superposition of the two corresponding linear mode shapes i and j . An expansion of Eq. (24) for two modes results in which indicates that the final coupled quadratic mode shape component is [13] (29) Here, g ij is the unadjusted quadratic mode shape component resulting from applying Eqs. (29)-(33) on the superposed artificial mode shape.

Impact on the frequency domain flutter equation
The quadratic mode shape components add a generalized aerodynamic forces term to Eq. (5) which is linearly dependent on the generalized coordinate [12,19]. Each element of the generalized aerodynamic forces matrix may be expressed as where f 1 j is the transfer function vector of unsteady aerodynamic forces induced by harmonic oscillation of mode j and f 0 the steady aerodynamic forces vector. The second term in Eq. (37) may be added to the generalized stiffness matrix in Eq. (5) as a geometric stiffness term K g (M ∞ ) leading to with

Simulation models
The configuration studied in this work is a generic T-tail configuration according to [5] (Fig. 1). The tail design features untapered as well as unswept vertical and horizontal tail planes. The simplicity of the model facilitates a fast method development and provides an opportunity for a comparison with literature data. Furthermore, the model can easily be modified for parametric studies.

Structural model
The structural model properties are chosen according to [5] and listed in Table 1. The finite element model, illustrated in Fig. 2, is set up as a beam model with massless flexible beams at 25 % chord of the VTP and the HTP. The mass properties are defined by concentrated mass elements with an offset of 10 % of the chord length aft of the elastic axis. The leading edges and trailing edges are realized by structural grid points rigidly connected to the flexible beams. These grid points serve the purpose of an improved interpolation of the structural mode shapes on the CFD surface mesh. At the VTP root, the finite element node of the flexible beam is constrained in all six degrees of freedom.
With the properties as listed in Table 1, the structural free vibrational characteristics feature a VTP out-of-plane bending eigenmode at 2.85 Hz (Fig. 3a) and a VTP torsion eigenmode at 5.29 Hz (Fig. 3b). The first structural mode results in an HTP roll motion with low amount of yaw while the second mode yields an HTP inplane motion with low amount of roll.

CFD model
An unstructured computational mesh for nonlinear, inviscid flow simulation is used for this study. The farfield covers 50 chord lengths in front, left, right, below, and above the configuration and 150 chord lengths aft of it. The surface mesh is shown in Fig. 4.
The present work uses a variation of the HTP incidence angle for studying the effect of steady aerodynamic forces on T-tail flutter in combination with the extended modal approach. Symmetric meshes w.r.t. the xz-plane are desired, for what reason a surface mesh of half of the T-tail with the xz-plane as symmetry plane is generated first. A rotation of the HTP's surface mesh followed by a volume mesh deformation using radial basis functions and a mesh repair step for badly shaped volume cells is used for the generation of different CFD meshes for each HTP incidence angle. A mirroring step then yields the final CFD meshes for the study.
A convergence study of the CFD mesh is carried out for an angle of incidence of 0.0 • at Mach 0.4 and generalized aerodynamic forces as target values. Three meshes with varying mesh densities and the properties as listed in Table 2 are considered. Fig. 5 displays the real and imaginary GAF matrix entries against the reduced frequency computed for the three mesh densities. As there is no distinct impact on the GAF matrix values, the computational mesh with medium density is selected.

Aeroelastic model
Deformations of the structure are interpolated on the CFD model using a moving least-squares (MLS) approach. The CFD surface meshes of the interpolated mode shapes are displayed in Fig. 6 (orange surface color) against the undeformed geometry (gray surface color).  Out-of-plane bending stiffness, EI 1 / N m 2 1.0E+07 1.0E+10 Inplane bending stiffness, EI 2 / N m 2 1.0E+10 1.0E+10

Quadratic mode shape components
The quadratic mode shape components related to the VTP bending and torsion modes are calculated using the nonlinear finite element approach described in Sect. 2.4.1 and the linear finite element approach described in Sect. 2.4.2. Fig. 7a-c illustrate these components calculated with the linear approach and already interpolated on the CFD surface mesh. The only notable difference between the linear and nonlinear finite element approach lies in the coupled quadratic mode shape component, for what reason Fig. 7d is inserted. The uncoupled quadratic mode shape component of the first VTP out-of-plane bending mode (Fig. 7a) is both, a shortening of the HTP and the VTP. The uncoupled quadratic mode shape component of the first VTP torsion mode is limited to a shortening of the HTP (Fig. 7b). The deformation of the coupled quadratic mode shape (Fig. 7c, additionally scaled) is a comparatively small shortening of the VTP. The nonlinear finite element approach yields an additional translational component in longitudinal direction along with a change in HTP incidence angle (Fig. 7d, additionally

Generalized stiffness
As described in Sect. 2.4.3, the additional generalized aerodynamic forces due to quadratic mode shape components are linearly dependent on the generalized coordinate and, hence, equivalent to stiffness terms. Fig. 8 illustrates this effect by plotting the change of the individual elements of the generalized stiffness matrix K arising from the extended modal approach against the steady HTP incidence angle (i.e., the steady aerodynamic HTP load). This figure is also suitable for comparing the two methods described in Sect. 2.4.1 ("NonlinearFE") and 2.4.2 ("LinearFE"). It can be observed that both methods agree well with each other in terms of the impact on the generalized stiffness matrix as well as the global trend w.r.t. the steady HTP incidence angle.
Regarding the modal stiffness of the first VTP bending mode ( K 11 , upper left plot in Fig. 8), one can recognize a reduction of the modal stiffness at negative HTP incidence angles (downforce at the HTP), while a positive HTP incidence angle (upforce at the HTP) increases this value. The absolute change in generalized stiffness is assigned to the first vertical axis, while the relative change w.r.t. the initial generalized stiffness without quadratic mode shape components is assigned to the second vertical axis. A relative change in generalized stiffness of roughly ±10 % is observed.
An order of magnitude smaller is the relative change in stiffness related of the VTP torsion mode ( K 22 , lower right plot in Fig. 8) with values up to 1 %. However, both, a negative incidence angle and a positive one, increase the stiffness.
The additional symmetric off-diagonal stiffness values ( K 12 and K 21 , upper right and lower left plots in Fig. 8) show a stiffness increase with increasing HTP incidence angle, similar to the additional stiffness due to the first VTP bending mode's quadratic mode shape component. As there are no off-diagonal generalized stiffness values from the linear modal analysis, the relative change is not depicted.

Flutter velocities
The flutter velocities of the generic T-tail configuration are computed using unsteady aerodynamic forces obtained from the LFD CFD solver TAU-LFD (Sect. 2.3). A non-matched approach is chosen with the aerodynamic forces being computed at a Mach number of 0.4 and atmospheric conditions acc. to the International Standard Atmosphere at mean sea level. In Fig. 9, the flutter velocities are plotted against the HTP incidence angle, where a negative incidence angle indicates again a downforce and a positive one an upforce at the HTP. A decreasing flutter velocity is observed with increasing HTP incidence angle up to approximately 1.0 • ("LFD", blue solid line). A further increase in HTP incidence angle then leads to a sudden increase in flutter velocity as well.
The uncoupled quadratic mode shape components ("LFD with gii (...)", yellow solid and red dashed lines) result in larger flutter velocities at negative incidence angles with deviations up to 2 %. A marginal decrease at positive incidence angles by less than 1 % is observed.
Adding the off-diagonal stiffness terms related to the coupled quadratic mode shape components ("LFD with gij (...)", green solid and purple dashed lines) results in a diminishing effect of the quadratic mode shape components for negative incidence angles w.r.t. the approach with only the linear mode shapes ("LFD", blue solid line). However, the effect of the off-diagonal stiffness terms becomes more pronounced at positive incidence angles with a considerable increase in flutter velocity. Fig. 8 Effect of quadratic mode shape components on generalized stiffness values of the generic T-tail configuration with variation of steady HTP incidence angle Fig. 9 Effect of quadratic mode shape components on flutter velocities of the generic T-tail configuration with variation of steady HTP incidence angle at Mach 0.4

Quadratic mode shape components
In simple terms, the first VTP bending mode can be regarded as a rigid body roll mode w.r.t. the VTP root in the yz-plane (Fig. 10a). The deformation of an arbitrary point r on the T-tail can then be described as Approximating the cosine and sine terms in the applied rotation matrix as Taylor polynomials [33] and retaining only the linear and quadratic terms leads to Applying this formalism to the T-tail model yields the deformations shown in Fig. 10. The linear deformation (first term in Eq. (44), Fig. 10a) is similar to the mode shape shown in Fig. 6a while the quadratic contribution (second term in Eq. (44), Fig. 10b) resembles the quadratic mode shape component shown in Fig. 7a.
In a similar manner, the first VTP torsion mode may be regarded as a rigid body yaw mode w.r.t. the VTP root in the xy-plane. The linearized deformation and the quadratic component are shown in Fig. 11. They are comparable to Fig. 6b and 7b , respectively. An additional chordwise shortening can be observed in Fig. 11b, but as the structural model does not feature flexible line elements in chordwise direction (see Sect. 3.1 and Fig. 2), this cannot be seen in Fig. 7b.
The analytical description of the coupled quadratic mode shape component is more extensive. Assuming the T-tail roll mode has a small yaw component and the yaw mode has a small roll component, the rotation matrices can be expressed as where () 1 indicates the predominant roll mode and () 2 the predominant yaw mode, while is the roll and the yaw angle. The coupled quadratic mode shape component is a superposition of the two rotation matrices. By retaining only the coupled quadratic terms in the superposed rotation matrix, the resulting quadratic deformation component becomes It can be seen that a coupling between a pure roll and a pure yaw mode ( 1 , 2 = 0 ) yields no quadratic deformation component. With non-pure roll and yaw modes, the quadratic deformation component is limited to a longitudinal (45)  deformation, which is proportional to the vertical coordinate, and a vertical deformation, which is proportional to the longitudinal coordinate. The linearized coupled rotation and its quadratic deformation component are shown in Fig. 12. Note the similarity between Fig. 12b and 7d.
Regarding the VTP bending and torsion modes of the generic T-tail unit, the coupling results from the non-overlapping mass and elastic structural axes. This is expected to be more pronounced for common T-tail units with additional VTP and HTP sweep and taper.

Generalized stiffness
Again, one can start with approximating the first VTP outof-plane bending mode with a rigid body roll mode around the VTP root, which has a rotational stiffness k . As shown in [4], an external follower force L 0 , acting in direction of the VTP, and a linear displacement approximation leads to the equation of motion where J x is the T-tail unit's mass moment of inertia w.r.t. the longitudinal axis and h the VTP span. Adding the quadratic deformation component compensates the spurious stiffness term in Eq. (48) (second term in brackets). Thus, the quadratic deformation component and the steady aerodynamic load at the HTP increase the stiffness proportionally to the span of the VTP and the steady load, similar to the results shown in Fig. 8. It should be noted here that the lack of symmetry in the change of stiffness observed for the generic T-tail configuration is related to the non-symmetric steady lift forces around a zero HTP incidence angle, as the zero-lift incidence angle is positive.
The increasing stiffness of the VTP torsion mode can be explained similarly, but with the lateral forces acting at the HTP tips and a dependence on the HTP semi span instead of the VTP span. As these forces are symmetric in sign and (48) J ẍ− (k − L 0 h) = 0 magnitude w.r.t. the zero-lift incidence angle, the resulting impact on the stiffness is also symmetric.
As shown in Sect. 4.1, the coupled quadratic mode shape component is partly a shortening of the VTP. Hence, its effect on the stiffness of the two linear mode shapes is comparable to that of the quadratic mode shape component due to VTP bending. Since the nonlinear FE approach also yields a longitudinal deformation component, the stiffness difference between the two approaches at higher HTP incidence angles may be related to steady drag forces. However, a conclusive evaluation is still to be carried out.

Flutter velocities
Regarding the influence of the uncoupled quadratic mode shape components on T-tail flutter, the decreasing VTP bending stiffness at negative HTP incidence angles due to the influence of the related uncoupled quadratic mode shape component (cf. Fig. 8) leads to a larger frequency ratio between VTP bending and torsion. This is amplified by a simultaneously increasing stiffness of the VTP torsion mode. Fig. 13 depicts the modal frequency ratio between VTP torsion and bending with uncoupled quadratic mode shape components and variation of the steady HTP incidence angle. An increasing value on the vertical axis indicates a larger separation of the frequencies. Since the unsteady aerodynamic forces are not modified by the extended modal approach, the flutter velocity increases with increasing frequency ratio as well. Vice versa, the flutter velocity decreases for positive HTP incidence angles due to the lower frequency ratio.
The off-diagonal stiffness terms originating from the coupled quadratic mode shape components introduce a static mode coupling. A softer coupling at negative HTP incidence angles results in a flutter velocity reduction, while a stiffer coupling increases the flutter velocity. A further study of this effect is subject of future activities.   13 Modal frequency ratio between VTP torsion and bending mode with uncoupled quadratic mode shape components and variation of steady HTP incidence angle

Conclusion and outlook
The effect of quadratic mode shape components on the flutter velocities of a generic T-tail configuration is shown using an extended modal approach and aerodynamic forces from LFD CFD. Two methods based on linear and nonlinear finite element analysis, respectively, are implemented to obtain the quadratic mode shape components of the two fundamental linear mode shapes of the T-tail, namely VTP out-of-plane bending and torsion. The g-method is employed to compute flutter velocities at various HTP incidence angles, i.e., steady HTP loadings.
At negative HTP incidence angles, the uncoupled quadratic mode shape components reduce the stiffness of the VTP bending and increase the stiffness of the VTP torsion mode. This results in an increase in flutter velocity. At positive HTP incidence angles, both T-tail mode shapes become stiffer with a more dominant stiffening of the VTP bending mode, which leads to a marginal decrease in flutter velocity.
The coupled quadratic mode shape component counteracts the effect of the uncoupled quadratic mode shape components for negative HTP incidence angles, leading to almost the same flutter velocities as if no quadratic mode shape components had been included. At positive incidence angles, however, the uncoupled quadratic mode shape component increases the flutter velocities substantially.
Subsequent studies will explore the nature of the differences between the linear and the nonlinear finite element approach in calculating coupled quadratic mode shape components. In addition, the impact of the off-diagonal stiffness terms on the flutter velocities in dependence on the static aerodynamic reference will be addressed. Moreover, the extended modal approach requires a verification by means of modal as well as fully coupled time domain simulations.