Abstract
The vertical tail buffet induced by the vortex breakdown flow is numerically investigated. The unsteady flow is calculated by solving the RANS equations. The structural dynamic equations are decoupled in the modal coordinates. The radial basis functions (RBFs) are employed to generate the deformation mesh. The buffet response of the flexible tail is predicted by coupling the three sets of equations. The results show that the presence of asymmetry flow on the inner and outer surface of the tail forced the structural deflection offsetting the outboard. The frequency of the 2nd bending mode of the tail structure meets the peak frequency of the pressure fluctuation upon the tail surface, and the resonance phenomenon was observed. Therefore, the 2nd bending responses govern the flow field surrounding the vertical tail. Finally, the displacement of the vertical tail is small, while the acceleration with a large quantitation forces the vertical tail undergoing severe addition inertial loads.
1. Introduction
LARGE swept delta wing combined with vertical tail configuration is widely used by modern flight vehicles. This kind of aerodynamic configuration is designed to reduce drag and improve the manipulated performance at a large attacking angle. However, the vortices emanating from the leading edge of the delta wing might breakdown before reaching the vertical tail at a certain attacking angle. In such cases, the vertical tails suffer from the wideband frequency excitation of the highly turbulent flow, which would produce a severe buffet on the tails. Wentz [1] conducted a wind tunnel experiment to investigate F/A-18 tail buffet problem, and the results revealed that the buffet is a phenomenon of resonance induced by the excitation of vortex breakdown flow, while the frequency bandwidth of the excitation covers the natural frequency of tails. The characteristics of buffeting are different from the traditional flutter or limit cycle oscillation (LCO). The buffet displacement with slight amplitude produces large acceleration responses that make the tail suffer from severe inertia force loads and fatigue failure.
Numerical studies for the aeroelastic responses in the vortical flow conducted before are mainly focused on the aerodynamic load, flutter, or LCO responses for lightweight materials. However, studies on the time/frequency characteristic of buffet responses are relatively inadequate. Findlay [2] simulated leading-edge vortex breakdown flow nearby vertical tails by solving thin-layer Navier–Stokes equations. It was found that the computed pressure distributions on the surface of the tails agree well with the experimental measurements. Kang [3] calculated the vertical tail buffeting loads at large attacking angle and obtained the buffeting characteristics. Kandil [4–8] investigated the bending-torsion interaction effect on the tail’s buffet, as well as the effect of freestream Mach number on the tail’s buffet. The results revealed that the deflections and loads of the coupled bending-torsion response case are substantially lower than those of the uncoupled bending-torsion response case, and the buffet loads in the transonic flow are lower than those in the subsonic flow. Guillaume [9] simulated the F/A-18 tail buffet in the transonic flow and obtained the displacement responses of the tails. Ahmed Elmekawy [10] conducted numerical simulations for F/A-18 tail buffet with four turbulence models and compared the computational results with the experiment data. It was found that the buffet problem could be simulated as a two-way fluid-structure interaction, and by using the nonlinear eddy viscosity model, since these give better results when compared to the other considered models. Attar [11] conducted numerical simulations for the full-span delta wing buffet in the vortex breakdown flow. Attar believes the lift enhancement phenomena are due to the reorganization of the flow and reformation of a leading-edge vortex structure.
The present study conducts a numerical simulation to investigate the tail buffet characteristics. Roe [12] proposed a finite different scheme for obtaining numerical solutions to Riemann problems. Results show that this fast scheme is very accurate and might be used to reduce the number of iterations. This method was employed in this paper to discretize the Euler flux terms of N-S equations. Yoon [13] developed a vectorizable and unconditionally stable scheme to solving the Euler and N-S equations. Results show that this method is efficient and robust. This method was used for the time integration of the N-S equations in this paper. Rodden [14] demonstrated a fast interpolation scheme to fitting the payload data of an infinite plane. This method was used here to perform the transformation of structural deformation from the CSD grid points to the CFD grid points. Jakobsson [15] proposed an interpolation scheme to propagate the deformations from the boundaries to the interior of the CFD mesh. Both structured and unstructured meshes can be treated by this method with low computational costs. This method was employed here to generate the deformation mesh. The unsteady Reynolds-averaged Navier–Stokes (RANS) equations with k-w SST turbulent model were used to model the vortical flow. The unsteady RANS solver was validated using the data from the Hummel’s 76-degree delta wing [16]. A structural model of the vertical tail was established by the finite element analysis. The flexible tail was positioned after the leading-edge vortex breakdown point to investigate the buffet characteristics. The deflection and acceleration responses are calculated by solving RANS equations coupled with structural dynamic equations.
2. Numerical Method
The numerical method for simulating the tail buffet was used to solve the aeroelastic equations, which consisted of the Navier–Stokes equations and structural dynamic equations. These equations coupled and advanced in time. Figure 1 shows the time stepping method.
Details of the methods are described as follows.
2.1. RANS Equations
The unsteady compressible Navier–Stokes (N-S) equations in the conservative form can be written aswhere t V represents the computational domain and V represents the computational boundary. The vector of conservative variables , J represents the Jacobian determinant, is density, u, , represents the velocity magnitude in three coordinate directions, e represents the total energy per unit mass, and t is time. All the variables mentioned above are dimensionless. A semidiscrete finite-volume formulation is used to discretize the N-S equations. The viscous flux terms are discretized by the central difference scheme, and the inviscid flux vector is discretized by FDS method of Roe [12], where the subscript variable i = 1, 2, 3, represents three coordinate directions. The LU-SGS implicit scheme [13] is used for time integration. The turbulence model was the k-ω SST model, which solves a two-equation turbulence model for the kinetic turbulent energy k and specific dissipation rate ω:where Gk represents the generation of turbulence kinetic energy due to mean velocity gradient. Gω represents the generation of ω. Yk and Yω represent the dissipation of k and ω due to turbulence. Dω represents the cross-diffusion term. Sk and Sω are user-defined source terms. σk and σω are the Prandtl numbers for k and ω. μt is the turbulence viscosity:where S is the strain rate magnitude. The modal constant a1 = 0.31, Ret = ρ·k/(μ·ω), and the blending function F is given bywhere y is the distance to the closest surface. Note that the constant coefficient α in the k − ω SST model is defined as α = α1·F + α2·(1 − F); the model constant α1 and α2 are 5/9 and 0.44, respectively.
2.2. Structural Dynamic Equations
The structural dynamic equations in the form of finite element analysis can be written as [11]where [M] is the mass matrix, [C] is the damping matrix, and [K] is the stiffness matrix. {U} is the displacement vector of nodal degrees of freedom and {F} is the aerodynamic force vector, both can be expressed aswhere N is the total number of node points of the structural model; ui and fi are vectors with three components in x, y, z directions:
In a modal approach, the modal decomposition of the structure motion can be expressed as the following:where [Λ] is the eigenvalue matrix, [Λ] = diag [λ1, λ2, … λ3N], λj = ωj2, where ωj is the natural frequency of jth mode, and the mode shape matrix [ϕ] = [ϕ1, ϕ2, …, ϕ3N]. Equation (8) can be solved by using a finite element solver (e.g., MSC.NASTRAN) to obtain [ϕ] and [Λ]. Thus, the displacement of the structure can be calculated as follows:
The energy of high order modes is relatively small. In this paper, only the first four modes are used in equation (9). Thus, the vector [a] = [a1, a2, a3, a4]T. Substituting equation (9) to equation (5) yields
Multiplying equation (10) by [ϕ]T yieldswhere , , , and .
Let ; equation (11) can be transformed in a state form as follows:where I represents the unit matrix. Equation (12) can be solved by the four-stage Runge–Kutta method. The infinite plane spline (IPS) method [14] was adopted as CFD-CSD interface data transformation method.
2.3. Mesh Deformation Formulations
The radial basis function (RBF) [15] was used to transform the CFD mesh. The interpolation function is defined aswhere F (r) denotes the interpolation function; N refers to the number of RBFs. In this paper, N is equal to the total number of the grid points on the surface. denotes the RBF, refers to the weighting factor of the function , r refers to the position vector of grid point, ri refers to the central point position vector of the RBF, and d refers to the influencing radius of the interpolation function. The Wendland’s C2 function was employed as :
The weighting factor can be obtained by solving the following equations:where is a vector of the displacement magnitude of each grid point on the surface, the subscript i = 1, 2, 3 denotes the x, y, z direction, respectively, the weighting factor vector W = [, ,…, ]T, and the coefficient matrix Φ,
3. Results and Discussion
The computational configuration is a half-delta wing with a vertical tail, as shown in Figure 2. The x, y, z coordinate, the wing chord c, and spanwise length of the wing L are described in Figure 2. More details of the model geometry and the structural parameters are described in Section 3.1. The wind tunnel experiments of the rigid delta wing are employed to verify the computation code. The buffet responses of the flexible tail as well as the instantaneous flowfields are demonstrated in Section 3.2.
3.1. Model, Grid, and Verification Test
The geometry of the delta wing is shown in Figure 3. The aspect ratio of the wing A = L2/S = 1, where S represents the wing area. The leading-edge swept angle of the wing is 75.96°. The root chord length of the wing is 1 m. The maximum thickness of the wing is 0.021 m. The vertical tail is located at the 1.06 chord length behind the apex of the wing, with a leading-edge sweep of 62.5°. The structure model of the tail was a single aluminum spar hidden inside a balsa-wood covering, as shown in Figure 4. The aluminum spar has a constant thickness of 0.0017 m. The root chord of the aluminum spar is 0.04 m, and the tip chord is 0.012 m. The aluminum spar was constructed from 6061-T6 alloy with a density of ρ = 2693 kg/m3, a moduli of elasticity of E = 6.896 × 1010 N/m, and a moduli of rigidity of G = 2.5925 × 1010 N/m. The root chord of the balsa-wood covering was 0.232 m, the tip chord was 0.054 m, and the spanwise length was 0.204 m. The thickness of the balsa-wood gradually decreases from 0.02 m at the tail root to 0.01 m at the tail midspan. From the midspan to the tail tip, the thickness was constant. The tail cross-section has a semidiamond shape with a bevel angle of 20°. The balsa-wood density was ρ = 179.7 kg/m3, the moduli of elasticity was E = 6.896 × 1010 N/m, and the moduli of rigidity was G = 2.5925 × 1010 N/m. The fundamental frequencies of the first and second bending mode of the tail are 33.941 Hz and 152 Hz, respectively, while the torsional frequency was 241.09 Hz, as shown in Figure 5. Figure 6 shows the CFD-CSD interface interpolation effect for the mode shapes of the vertical tail.
(a)
(b)
(c)
(a)
(b)
(c)
An O-H grid of 87 × 107 × 51 grids that points in the wrap-around, streamwise, and normal directions, respectively, was used for the solution. Figure 7 compares the original mesh nearby the delta wing/tail surface with the mesh deformation. The flow conditions were Rec = 1.25 × 106, the attacking angle was α = 40°, and the freestream Mach number was Ma = 0.3. The unsteady time step was chosen as 10−5 s to capture the high frequency vortex breakdown flow.
(a)
(b)
The experiment [16] for the vortex flow on the delta wing was employed to validate the computational capability of the CFD procedure. The freestream Mach number was 0.3, Rec = 2 × 106, and the attacking angle was α = 20.5°. The wing surface pressure measurement positions were x/c = 0.3, 0.5, 0.7, and 0.9. Figure 8 compares the computed surface pressure with the experimental measurements. It was shown that the computed solutions agree well with the experimental measurements. Figure 9 compares the calculated surface flow patterns with the experiments. It was shown that the main and the secondary separation lines were both captured by the numerical procedure. The results above show that the computational method presented in the present study is reliable.
(a)
(b)
3.2. Vertical Tail Buffet Results
The computational instantaneous streamlines at 0.2 s and 0.25 s are shown in Figure 10, and the total pressure slices and tip deflection of the vertical tail at the position of x/c = 1.3 were also extracted. It is found that the leading-edge vortices proceed to the outward side of the vertical tail, resulting in the tail deflection towards the outward side. The buffet load coefficients ΔCP nearby the tip of the tail (at the position of 90% spanwise length of the tail, 50% local root length of the tail) and the root of the tail (at the position of 30% spanwise length of the tail, 50% local root length of the tail) were both extracted in Figure 11. It was found that the mean value and the fluctuation magnitude of the buffet loads at the tip of the tail were both larger than those at the root of the tail, and the frequency bandwidth of the fluctuation loads was from 100 Hz to 200 Hz. In addition, the peak frequency of the buffet loads at both the tip and root of the tail appears at 153 Hz, which is almost the same as the fundamental frequency of the second bending mode of the tail, demonstrating the tail structure resonance under the fluctuation load excitation of the vortex breakdown flow.
(a)
(b)
(c)
(d)
Figure 12 compares the displacement responses with the acceleration responses on the midpoint of the leading-edge of the tail. It was observed that the amplitude of the displacement responses was small (magnitude: 10−4 m). However, the amplitude of the acceleration responses was 5∼6 orders of magnitude greater than the displacement responses (magnitude: 103 m/s2). Furthermore, the acceleration amplitude clearly diverges with time. The PSD curves in Figure 12 show that the peak frequency of both displacement and acceleration responses was 151.23 Hz, which meets the fundamental frequency of the second bending mode. The results above conclusively show that the buffet deflections with slight amplitude produce acceleration responses with a large amplitude, causing severe additional internal loads on the structure, and the resonance between the vortex breakdown flow and the second bending mode dominates the characteristics of the acceleration responses.
4. Conclusion
A computation method for solving RANS equations coupled with structural dynamic equations was presented in the present study to simulate the vertical tail buffet responses in the vortex breakdown flow. A validation computation of the vortex flow on Hummel’s 76° delta wing was conducted to verify the numerical method. A delta wing/vertical tail configuration was modeled and simulated at a large angle of attack to obtain the buffet response of the tail.
It was observed that the vortex breakdown flow on the inner and outer surface of the tail was asymmetry, and the vortices were located at the outer side of the tail, which makes the tail deflect to the outer side. The peak frequency of the buffet load was almost the same as the fundamental frequency of the second bending mode of the tail. The response amplitude of the second modal displacement and acceleration diverge with time. Meanwhile, the amplitudes of the first and third modal responses were irregular oscillations with time. The frequency of the modal displacement and acceleration responses was the fundamental frequency of each structural mode. In addition, the deflection amplitude of the tail is small (magnitude: 10−4 m). However, the acceleration amplitude was much greater (magnitude: 103 m/s2).
The results above conclusively show that the resonance between the vortex breakdown flow and the second bending mode dominates the characteristics of the acceleration responses, and the buffet deflections with slight amplitude produce acceleration responses with large amplitude, causing severe additional internal loads on the structure.
Data Availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the Shanxi Foundation Research Projects for Application (201801D221234) and Science and Technology Projects of Shanxi Province Guided by Central Administration of China (YDZX20201400001519).