Simulation on Buffet Response and Mitigation of Variant-Tailed Aircraft in Maneuver State

Abstract

This study proposes a computational fluid dynamics and computational structure dynamics (CFD/CSD) coupled method for calculating the buffet response of a variant tail wing. The large-scale separated flow in the buffet is simulated by the detached vortex approach, vibration deformation of the tail wing is solved by the dynamic mesh generation technique, and structural modeling is based on the mode method. The aerodynamic elastic coupling is calculated through the cyclic iteration of aerodynamics and the structural solution in the time domain. We verify the correctness of the proposed method through a typical delta wing calculation case, further simulate the buffet response of a variant tail wing in maneuver state, and finally realize buffet mitigation using an active excitation method. Overall, this study can provide an important reference for the design of variant-tailed aircraft.

1. Introduction

With the increasing demand for aircraft to perform complex tasks and adapt to different flight environments, variant technologies have received widespread attention [1]. Among all variant solutions, the folding wing or tail wing combines lift drag characteristics and stealth characteristics while also solving lateral stability and control problems, making it a highly promising variant technology. However, variant aircraft typically have low structural stiffness, complex mechanisms, and material compositions, and the significant dynamic changes in mass distribution, stiffness distribution, and aerodynamic load distribution during the variant process have brought about many new aeroelastic problems.

In recent years, many scholars have conducted exploratory research on the aeroelasticity of variant aircraft. John and Robert [2] integrated finite element modeling and analysis, an aerodynamic program based on the vortex lattice method, and the SIMULINK control module to conduct joint simulation on flexible Z-shaped wings. They simulated the dynamic response of folding wings under different wind tunnel test conditions and determined the influence of flight changes and deformation states on structural aeroelasticity. Reich et al. [3] used the IAMMS method to establish rigid and flexible full aircraft models of foldable wing aircraft and analyzed the changes in flight parameters such as the attack angle, rudder deviation angle, and velocity. Wilson et al. [4] and Castrichini et al. [5] studied the aeroelastic characteristics of a foldable wing tip aircraft, mainly exploring the effects of hinge orientation, hinge stiffness, hinge damping, and wing tip weight on the static and dynamic response, structural load, and structural weight of the aircraft. Lee and Chen established an approximate analytical nonlinear equation for [6] folding wing hinge structures and simulated the nonlinear aeroelastic characteristics of folding wing structures with free hinges. Tang et al. [7] and Wang et al. [8] studied the structural dynamics characteristics of foldable wing aircraft based on the foldable wing flat plate model to develop a deeper understanding of aerodynamic elastic stability and its potential effects; Zhao et al. [9] and Huang et al. [10] proposed a parameterized aeroelastic modeling method for folding wings using substructure synthesis and the dipole grid method. This method can achieve rapid buffet analysis of folding wings under different configurations. The analysis results show that the aeroelastic characteristics of wings under different folding angles will undergo significant changes.

The above research mainly focuses on the aeroelastic stability under fixed configurations and the dynamic response during variant processes, with little attention paid to the issue of buffet. Taking the variant tail wing aircraft that this article focuses on as an example, when in a high-speed heading control and combat state, the aircraft needs to maneuver. At this time, when the angle of attack is large, the detached vortices generated by the leading edge of the wing will rupture in front of the tail wing, forming a highly turbulent and rotating unsteady tail vortex flow, causing the tail wing to be immersed in highly turbulent vortices. The pulsating pressure caused by turbulence will stimulate the tail wing to oscillate. Due to the high excitation energy of the rupture vortex generated by a high angle of attack maneuvering flight, if the frequency band of this vortex covers the natural frequency of one or several modes of the tail wing, the tail wing will suffer from very serious buffet problems.

Currently, there are two main types of simulation methods for vibration response. One type is the calculation method, which uses the buffet pressure of a rigid wind tunnel model as the excitation spectrum and then applies frequency domain unsteady aerodynamic calculation methods such as the dipole grid method to account for the unsteady aerodynamic forces induced by wing vibration. Then, based on the random vibration spectrum calculation, the buffet response of the structure is estimated. Representative works include Lee [11] and Pototzky [12]. The advantage of this method is that the model and algorithm are simple, but there are significant limitations in accuracy and closeness to reality. Another type is the CFD-based buffet response calculation method, which can be divided into two types: unidirectional coupling and bidirectional coupling. The former believes that the separated flow/shock wave oscillation as an excitation is independent of structural motion, and the aeroelastic feedback effect is weak. Therefore, a decoupled two-step method is adopted for research. The unidirectional coupling method is still the main method used for studying the vibration response in engineering. The bidirectional coupling method (also known as CFD/CSD coupling method) achieves more realistic fluid–structure coupling simulation by iterating unsteady flow and structural dynamics calculations in the time domain. At present, this method has been widely developed and applied in the study of aerodynamic problems [13,14] such as buffet [15,16], but faces challenges such as high accuracy requirements for separated flow simulation and long response simulation time history requirements for buffet response problems. In limited research, Gardner et al. [17] coupled CFD based on Euler equations with finite element calculations based on Von Karman plate theory to predict the buffet response of a 60-degree sweep angle delta wing. Kandil et al. [18] coupled CFD based on RANS equations with structural modal calculations to study the structural response of triangular wing layouts with square vertical tails in different configurations. Sheta et al. [19] used a similar method to predict the vertical tail buffet response of F/A-18 aircraft.

Therefore, this study establishes a set of aerodynamic and elastic calculation methods for variant tail wing aircraft based on CFD/CSD coupling, focusing on unsteady flow DES (detached vortex) simulation, dynamic mesh generation, flow field structural field coupling, and other technologies to improve the fidelity and computational efficiency of buffet response prediction. The method was validated using a high angle of attack delta wing example, and buffet response simulation analysis was conducted on the maneuvering state configuration of a typical variant tail wing aircraft. The buffet response mitigation technology based on active excitation was explored.

2. Numerical Method

2.1. Unsteady Flow Simulation Method

For unsteady flow CFD simulation, the conservative form of the three-dimensional N-S equation in the Cartesian coordinate system using the arbitrary Euler–Lagrange description method is as follows:

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \left[\rho \right(\mathit{v}-{\mathit{v}}_g\left)\right]=0,$$

(1a)

$$\frac{\partial \rho \mathit{v}}{\partial t}+\nabla \cdot \left[\rho \mathit{v}\right(\mathit{v}-{\mathit{v}}_g)+p\mathit{I}]-\nabla \cdot \overline{\overline{\mathit{\tau }}}=0,$$

(1b)

$$\frac{\partial \rho E}{\partial t}+\nabla \cdot \left[\rho E\left(\mathit{v}-{\mathit{v}}_g\right)+pv\right]-\nabla \cdot \left(k\nabla T+\overline{\overline{\mathit{\tau }}}\cdot v\right)=0,$$

(1c)

In the above equations, ρ, p, E, k, and T are density, pressure, total energy per unit mass, thermal conductivity coefficient, and temperature. Respectively, they are fluid motion velocity vectors, grid motion velocity vectors, viscous stress tensors, and unit matrices.

Based on the strategy of spatiotemporal separation and discretization, a multi-block structured mesh finite volume method is used to integrate Equation (1), as follows:

$$\frac{\partial }{\partial t}{\displaystyle {\int }_{\Omega }W}d\Omega +{\displaystyle {\oint }_S(F_c-F_v)dS}=0,$$

(2)

Among them, S is the boundary interface of any control volume unit $\Omega$. W is a conserved quantity. ${\mathit{F}}_c$ and ${\mathit{F}}_v$ are the convective flux and viscous flux, respectively. Their specific forms can be found in [20].

In terms of spatial discretization, the convective flux calculation adopts the Roe scheme, which combines the high-precision characteristics of boundary layer capture and the high-resolution characteristics of shock wave capture. To obtain solutions with second-order accuracy or higher, a MUSCL format with Van Albada limiters was used to reconstruct the variables on both sides of the interface, ensuring the computational accuracy of the smooth region of the flow field and eliminating non-physical oscillations near the shock wave. To avoid interference between various physical quantities and reduce numerical oscillations, this study adopts a reconstruction method based on characteristic variables. The calculation of viscous flux adopts a second-order central difference scheme.

In terms of time discretization, due to the need for precise time advancement in the calculation of unsteady flow fields, this study adopts the dual time step method, which transforms solving unsteady problems into solving steady problems in the virtual time domain. This eliminates the limitation of physical time step size on stability conditions in the flow field and can still use the acceleration convergence measures in steady calculations. The time derivative term in Equation (2) is approximated using a three-point backward difference and the residual term is implicitly processed, resulting in

$$\frac{3{\Omega }^{(n+1)}{W}^{(n+1)}-4{\Omega }^{(n)}{W}^{(n)}+{\Omega }^{(n-1)}{W}^{(n-1)}}{2\Delta t}=-R^{(n+1)},$$

(3)

Among them, R is the residual vector; Δt is the physical time step size; and superscripts n − 1, n, and n + 1 represent the previous time, current time, and next time, respectively.

Introduce dimensionless time, as follows:

$${\Omega }^{(n+1)}\frac{d{W}^{\ast }}{d\tau }=-\left(R(W^{\ast })+\frac{3W^{\ast }}{2\Delta t}-Q^{\ast }\right),$$

(4)

Among them, W* and R(W*) are approximations of W⁽ⁿ⁺¹⁾ and R⁽ⁿ⁺¹⁾, respectively; Q* is a non-stationary source term.

After obtaining the solution to the above equation (Q* = 0 at that time), Equation (2) naturally satisfies and exists W⁽ⁿ⁺¹⁾ = W*. We determine the implicit time discretization of the above equation, and first write its first-order Euler implicit form (omit the asterisk representing the unsteady term) as follows:

$$\Omega \left(\frac{1}{\Delta \tau }+\frac{3}{2\Delta t}\right)\Delta W^{(n)}=-R^{(n+1)},$$

(5)

In the formula, ΔW⁽ⁿ⁾ = W⁽ⁿ⁺¹⁾ − W⁽ⁿ⁾ is the increment in the conserved quantity. Due to the unknown residual value at the next moment, we linearize it with respect to the current moment, resulting in

$$\left[\Omega \left(\frac{1}{\Delta \tau }+\frac{3}{2\Delta t}\right)I+\frac{\partial R}{\partial W}\right]\Delta W^{(n)}=-R^{(n)},$$

(6)

Solve the equation using LUSGS format to obtain ΔW⁽ⁿ⁾ to update the conserved variable ΔW⁽ⁿ⁾ at the next time step.

In terms of boundary conditions, the far field adopts non-reflective boundary conditions, the object surface adopts slip-free and adiabatic stagnation boundary conditions, and the symmetric surface adopts symmetric boundary conditions. To improve computational efficiency, parallel algorithms based on MPI are adopted in a multi-block structured grid system. The main idea is to divide the physical domain into several subdomains, and each subdomain’s computing task is separately assigned to a CPU core for calculation. During the calculation process, data communication between different cores is achieved through MPI, providing a unified interface function.

Unsteady flow in a buffet environment exhibits highly turbulent characteristics, and turbulence simulation is a key factor affecting the prediction of buffet response. Due to the significant limitations of the RANS method in simulating large-scale eddies and the enormous computational resources required by the LES (large eddy simulation) method, this study adopts a DES simulation method that combines the characteristics of both methods. The main idea is to solve the RANS equation within the boundary layer to simulate turbulent flow within the boundary layer while using the large eddy simulation method in the outer region.

The control equation of the DES turbulence model is based on the SA equation, as shown in Formula (7):

$$\frac{D\tilde{\nu }}{Dt}=C_{\mathrm{b}1}\tilde{\Omega }\tilde{\nu }-C_{\mathrm{w}1}{f}_{\mathrm{w}}{\left(\frac{\tilde{\nu }}{d-f_{\mathrm{d}}\max (0,d-C_{\mathrm{DES}}\Delta )}\right)}^2+\frac{1}{\sigma }\left\{\nabla \cdot \left[\left({\nu }_{\mathrm{L}}+\tilde{\nu }\right)\nabla \tilde{\nu }\right]+C_{\mathrm{b}2}{\left(\nabla \tilde{\nu }\right)}^2\right\},$$

(7)

In the formula, $\tilde{v}$ is the control variable of the turbulence equation, $v_L$ is the laminar viscosity coefficient, $\tilde{\Omega }$ is the local rotation rate, d is the closest distance to the object surface, Δ is the local grid scale, and the rest are empirical constants. When d « Δ, the above equation is equivalent to the RANS simulation; when d » Δ, it is equivalent to the Smagorinski large eddy simulation. The purpose of introducing the switch function $f_d$ is to eliminate the disadvantage of the calculation results in the RANS and LES boundary area being greatly affected by the computational grid [21]. An effective construction method is shown in Formula (8).

$$f_{\mathrm{d}}=1-\tanh (512{r_{\mathrm{d}}}^3),r_{\mathrm{d}}=\frac{{\nu }_{\mathrm{L}}+\tilde{\nu }}{\sqrt{U_{i,j}{U}_{i,j}}{\kappa }^2{d}^2},$$

(8)

Among them, $U_{i,j}$ represents the velocity gradient, κ is a constant, and $r_d$ is designed to approach 1 in the logarithmic region of the boundary layer, while gradually decreasing to 0 in the viscous outer region of the boundary layer, thereby achieving a reasonable transition between RANS and LES.

Discretize the turbulence control equation using the finite volume method for any control volume element Ω:

$$\frac{\partial }{\partial t}{\displaystyle {\int }_{\Omega }\tilde{\nu }d\Omega }+{\displaystyle {\oint }_{\partial \Omega }\left(F_{c,T}-F_{v,T}\right)dS}={\displaystyle {\int }_{\Omega }{Q}_{T}d\Omega },$$

(9)

In the above equation, $F_{c,T}$, $F_{v,T}$, and $Q_T$ represent the convective flux term, viscous flux term, and source term, respectively.

Discretize $F_{c,T}$ and $F_{v,T}$ using the first-order upwind scheme and central scheme, respectively, and the time discretization method is consistent with the conservation equation above. On the far field boundary, $\tilde{v}$ is taken as 3–5 times $v_L$; on the surface boundary, $\tilde{v}$ = 0 is taken, and the initial value is calculated as 3–5 times $v_L$. Because the real data in the far field are difficult to obtain, the value of $\tilde{v}$ is mainly based on practical experience and previous literature [21].

2.2. Dynamic Grid Generation Technology

The buffet of variant tail wing aircraft belongs to the dynamic boundary problem, and describing its complex deformation requires high requirements for dynamic mesh generation methods. Considering that a single moving grid method experiences difficulty in balancing computational efficiency and deformation ability, this study adopts a self-developed hybrid moving grid method, RBF-TFI [22], which combines the strong deformation ability and high interpolation accuracy of the RBF (Radial Basis Function) method, as well as the high computational efficiency and adaptability of the TFI (Transfinite Interpolation) method for large-scale computational grids. It can also achieve the automatic generation of dynamic grids for complex shapes and deformations.

For any directional component of the deformation, the RBF interpolation function f can be expressed as the weighted sum of a series of basis functions φ. In this study, the thin plate spline function φ(x) = x²log(|x|) is selected as the basis function:

$$f\left(x\right)={\displaystyle \sum _{n=1}^N{a}_n{\Vert x-x_n^{\ast }\Vert }_2^2\log \left({\Vert x-x_n^{\ast }\Vert }_2+\delta \right)},$$

(10)

The deformation on N control points is known, where x is the coordinate vector of any grid point, $x_n^{*}$ is the coordinate vector of the nth control point, $a_n$ is the coefficient to be solved, and δ is a small quantity.

The above equation can be rewritten as a matrix form for all control points. By using the Gaussian full principal element elimination method, $a_n$ can be obtained, and then the deformation of any grid point can be calculated.

Using the Boolean sum of interpolation functions, three-dimensional TFI can be written as

$$\Delta {\mathit{x}}_{i,j,k}=U+V+W-UV-VW-UW+UVW,$$

(11)

Among them, ${∆x}_{i,j,k}$ is the deformation amount of any grid point. All the definitions of single and composite variables in the above equation can be found in [23].

The main steps of the RBF-TFI method for multi-block structured grid systems are as follows:

(1)

Properly select feature points on the edge interface of the grid block (far-field, object surface, and symmetry plane) as the control points required for RBF interpolation;

(2)

Use the RBF interpolation method to solve the deformation of grid block edges based on the deformation of control points;

(3)

Use the two-dimensional TFI method to solve the deformation of the grid block surface based on the deformation of the grid block edges;

(4)

Use the three-dimensional TFI method to calculate the deformation of all internal points of the grid block based on its surface deformation.

The practice has shown that for the variant tail wing aircraft used in this article, the RBF-TFI method can automatically generate dynamic meshes in seconds, and the quality of the deformed mesh is comparable to the initial mesh.

2.3. Flow Field Structural Field Coupling

Assuming that the aircraft structure is linear, the concept of vibration mode is used to describe the structural deformation Δr, as follows:

$$\Delta r\left(x,y,z,t\right)={\displaystyle \sum _{i=1}^n{h}_i\left(x,y,z\right)q_i\left(t\right)},$$

(12)

In the above equation, n is the selected number of modes, $h_i$ is the i-th mode vector, and $q_i$ is the generalized coordinate of the i-th mode.

Using the Lagrangian equation, the structural motion equation can be described in matrix form, as follows:

$$M{q}_{tt}+G{q}_t+Kq=A,$$

(13)

Among them, M, G, and K represent generalized mass, generalized damping, and generalized stiffness matrices, respectively. The generalized aerodynamic force A is related to the unsteady distributed load on the surface of the aircraft, as follows:

$$A_i={\displaystyle ∯\left[p(x,y,z,t)+\tau (x,y,z,t)\right]h_i(x,y,z)dS},$$

(14)

Among them, p(x, y, z, t) and τ(x, y, z, t) represent the unsteady pressure and friction stress distribution on the surface of the aircraft, respectively, obtained through the above unsteady CFD calculation.

The structural vibration modes have orthogonality, and after mass normalization, M and K are both diagonal matrices that satisfy

$$K_{ii}={\omega }_i^2{M}_{ii},$$

(15a)

$${\omega }_i=2\pi f_i,$$

(15b)

Introduce state variable E to reduce the order of the structural motion equation, as follows:

$$E={\left[q_1,q_2,\cdot \cdot \cdot ,q_n,q_{t1},q_{t2},\cdot \cdot \cdot ,q_{tn}\right]}^T,$$

(16)

After substituting into Equation (13), the following linear system of equations is obtained:

$$E_t=\left[\begin{array}{cc}0& I\\ -M^{-1}K& -M^{-1}G\end{array}\right]E+\left[\begin{array}{c}0\\ M^{-1}\end{array}\right]A,$$

(17)

Among them, 0 is an n-order zero square matrix, and I is an n-order unit matrix. The second-order estimation correction method [24] is used to solve the above equation system, thereby achieving efficient coupling between aerodynamics and structure.

Due to the different physical properties of the fluid domain and the solid domain, the computational grids of the two domains usually do not coincide on the coupling boundary. Therefore, it is necessary to interpolate the displacement at the structural points to the aerodynamic points. This study adopts the Infinite Plate Spline (IPS) interpolation method [25], whose core idea is based on the principle of virtual work. By introducing the concept of load vector corresponding to the vibration mode vector, the differential control equation of the thin plate is solved to achieve displacement interpolation. Traditional aeroelastic research requires repeated implementation of interpolation from structural deformation to flow field and from flow field load to structure, and force interpolation requires additional conservation conditions compared to displacement interpolation, resulting in higher complexity. To this end, the concept of spatial mode is introduced, and the structural vibration mode is interpolated to the aerodynamic surface grid points. The dynamic grid generation method is used to extend the vibration mode of the aerodynamic surface grid points to the entire flow field spatial grid, thereby improving the efficiency of flow field structural field coupling.

Figure 1 shows the detailed calculation process of the buffet response of a variant tail wing aircraft.

3. Calculation Results and Analysis

Firstly, an example of a delta wing buffet was used to verify the accuracy of the calculation method proposed in this study. Then, the buffet response of a typical variant tail wing aircraft was simulated to determine its maneuvering state configuration. Finally, a buffet response mitigation technique based on active excitation was explored.

3.1. Verification Example of High-Angle-of-Attack Buffet Response of Delta Wing

The calculation example is a 50-degree sweep angle flat delta wing [21] (as shown in Figure 2), with the calculated state being the inflow velocity of 30 m/s, angle of attack of 20°, and Reynolds number of 6.2 × 10⁵ based on root chord length.

In this study, the flow field and structural response of the delta wing are assumed to be symmetric; therefore, a semi-model calculation is considered. The flow field calculation grid adopts an H-H topology structure, with a total of 6.4 million grids divided into 48 blocks for parallel calculation. Increasing the grid density of the vertex, leading edge, and trailing edge of the delta wing, and the height of the first layer of mesh near the object surface is about 1 × 10⁻⁵ m, resulting in y⁺ ≈ 1. The flow field calculation grid is shown in Figure 3.

The material of the triangular wing structure is aluminum, with a density of 2700 kg/m³, a Young’s modulus of 6.9 × 10¹⁰ Pa, and a Poisson’s ratio of 0.3. The finite element modeling adopts a four-node shell element, with a total of 851 elements. The root adopts fixed boundary conditions, and the structural field calculation mesh is shown in Figure 3b. The first five modes obtained from structural finite element analysis were selected for buffet response simulation. A comparison between the calculated natural frequencies and the results in [26] is shown in Table 1, and the two are in good agreement.

Figure 4 shows the first five mode shapes of the delta wing structure, and their shapes are consistent with those in the literature.

Figure 5 shows the calculated non-dimensional wing tip displacement response Δz/c and acceleration response Δ$\ddot{z}$, where c is the half-span length. As shown in the figure, after a brief transition phase, a stable structural vibration response appeared, indicating that the delta wing experienced a buffet phenomenon at this angle of attack.

To quantitatively compare the vibration response, a time-domain characteristic analysis was conducted on the vibration response time history, including the displacement mean, displacement amplitude, and root mean square of acceleration response. Table 2 presents a comparison of the predicted characteristics of this project with the literature [26] and experimental measurement results [21]. It can be found that the mean displacement response obtained in this article is smaller than the experimental value, while the amplitude matches well, which is consistent with the literature results. For the root mean square of the acceleration response, the predicted values in this study are relatively close to the experimental values.

As shown in Figure 6, the fast Fourier transform (FFT) plots of the displacement response and acceleration response were obtained. Regarding the displacement response, the maximum peak appears around 21 Hz, while the peak values are smaller near the second-order, third-order, and higher-order frequencies, indicating that the dominant displacement vibration form at this angle of attack is the first-order bending mode. Unlike the displacement response, the acceleration response has significant peaks near the natural frequencies of each mode, and the peak sizes are generally similar. Similar phenomena also occur in [26].

Figure 7 shows the calculated instantaneous flow field structure, which is visualized using vorticity cloud maps at typical chord positions. It can be found that vortex breaking occurs near the vertex of the delta wing, which then develops into a large-scale unsteady flow field region, while the shedding vortex structure near the leading edge approximates complete stall flow. This highly unsteady flow field generates significant pressure disturbances on the wing surface, thereby stimulating the structural dynamic response of the delta wing.

3.2. Simulation of Buffet Response of Variant-Tailed Aircraft

A small aspect ratio variant tail wing aircraft is taken as the research object, and a typical maneuvering state is selected for buffet response simulation research. The reason is that in this state, the aircraft is subjected to a large angle of attack of the airflow; at this time, the detached vortices of the front wing and their rupture may cause the tail wing buffet phenomenon. The total length of the aircraft is 18.92 m, with a span of 14.50 m. The wing sweep angle is 45°, and the tail wing is 45° upward. Calculate the state at Mach 0.5, with a speed and pressure of 17.7 Ka, and a variable angle of attack.

The calculation of unsteady flow field adopts H-H type body-fitted structured mesh, with a mesh size of about 3 million. The calculation mesh is shown in Figure 8.

Firstly, the unsteady flow field calculation at high angles of attack was carried out for the rigid model of the variant tail, which helps to understand the mechanism of buffet response and also provides a preliminary analysis of the critical angle of attack at which buffet occurs. Figure 9 shows the instantaneous flow field at different typical angles of attack, and the results show that as the angle of attack increases, the intensity of the vortex increases, and the position of the vortex core is closer to the inner side of the fuselage. At an angle of attack of 20°, the range of immersion of the tail wing in the vortex is maximum, which is likely to induce buffet, and the potential buffet response is maximum. When the angle of attack significantly decreases or increases, the intensity of the buffet response will decrease.

Next, buffet response calculation was carried out for the elastic model of the variant tail wing while considering that the fuselage and wing of the aircraft are rigid, while the tail wing is elastic. The finite element model of the tail wing structure is shown in Figure 10.

The modal method was used for structural modeling, and all calculations took the first six modes with natural frequencies of 27.6, 48.5, 196.9, 273.5, 327.7, and 483.7 Hz, respectively. The first two vibration modes are shown in Figure 11, which are the first bending mode and the rotating mode. The bending mode of the structure is shown in Figure 11a, which is an elastic mode, and the rotating mode is shown in Figure 11b, which is a rigid mode. Two single-point fixed constraints are applied to the tail root to simulate the connection between the tail and the fuselage.

Based on the CFD/CSD coupling method established above, Figure 12 shows the generalized displacement response curves at different angles of attack. Due to the small response of the third order and above, only the calculation results of the first and second modes are shown in the figure. The period from the start to 0.2 s is the transition calculation phase of the response.

In Figure 12, it can be found that the generalized displacement response at low angles of attack exhibits a stable and regular attenuation, indicating that there has not been any chattering phenomenon. The response at α = 10° is continuously attenuated, but the attenuation rate is small, and the frequency is about 46 Hz, close to the frequency of the rotating mode, indicating that there is a small damping phenomenon in the rotating mode, which may be related to the neglect of structural damping. At α = 15°, the first-order mode exhibits an irregular and non-decaying response and begins to exhibit a chattering phenomenon. As the angle of attack increases further, the chattering response first increases and then decreases, reaching a peak at α = 20°.

Figure 13 further illustrates the wing tip displacement response curves at different angles of attack, which provides a more intuitive view of the maximum vibration response intensity at α = 20°. For a clear comparison, the response calculation results in Figure 13 are presented from 0.2 s.

To further analyze the vibration response, a spectral analysis of the displacement response of the wing tip was conducted after stabilization at α = 20°, and Figure 14 shows the power spectral density results. It can be found that the maximum peak of the response occurs around 28 Hz, and at this angle of attack, the dominant vibration form is the first-order bending mode. The second two peaks appear around 195 Hz (third mode) and 327 Hz (fifth mode), respectively, while the peaks near other modes are very small.

3.3. Simulation of Buffet Response Reduction for Variant-Tailed Aircraft

Currently, measures to mitigate the vibration response can be divided into two categories: passive and active. The former mainly starts from an aerodynamic perspective, reducing the generation of vortices by changing the shape or guiding the forebody vortices to change direction to avoid the tail wing. The latter mainly starts from a structural perspective and suppresses the vibration response by introducing devices such as vibration dampers and piezoelectric actuators. This study explores the vibration reduction technology based on piezoelectric actuators and simulates it by active excitation force on the wing surface. Due to the dominant vibration form of the buffet response being the first-order bending mode, the excitation position is selected near the half-chord length of the wing tip, and the excitation direction is selected as the local surface normal.

According to the principle of mechanical vibration, this study designs an excitation force f = −bΔ$\dot{x}$ based on velocity feedback, where Δ$\dot{x}$ is the instantaneous normal velocity at the excitation point, and b is the velocity feedback coefficient, which determines the magnitude of the excitation force while considering three scenarios: b = 0.1$b_0$, b = $b_0$, and b = ${10b}_0$, where $b_0$ is the reference coefficient that is related to the flight status, structural parameters, etc.

Figure 15 shows the displacement response curves of the wing tip under different excitation forces. It can be found that after active excitation, the amplitude of the buffet response decreases, and the greater the excitation force, the more significant the slowdown in the buffet response.

The power spectral density distribution of the wing tip displacement response is shown in Figure 16, which provides a more intuitive understanding of the damping effect of the active excitation on the buffet response. When the excitation force coefficient is set to b = ${10b}_0$, there is no obvious response peak.

4. Conclusions

This paper presents a numerical simulation study on the buffet response of a variant tail wing aircraft using the CFD/CSD time-domain coupling method. Firstly, key technologies, such as unsteady flow DES simulation, RBF-TFI dynamic mesh generation method, and flow structure field coupling, were introduced. Secondly, the developed method was validated using a high angle of attack delta wing example; based on this, a buffet response simulation analysis was conducted on the maneuvering state of typical variant tail wing aircraft. Finally, a simulation of active damping of the vibration response based on the application of excitation methods was also carried out. The research results of this study provide theoretical guidance and technical support for the design of future variant aircraft.