Fully Coupled Aeroelastic Analyses of Wing Flutter towards Application to Complex Aircraft Configurations

In previous studies (Yuan et al. 2013, 2016), a linear two-degree-of-freedom (two-DOF) structural model has been developed for the laminar-separation flutter problem at transitional Reynolds numbers and for transonic flows, respectively. The 2D two-DOF CSD module was transplanted and incorporated with the FLOWer CFD solver in this study. The governing equations for the 2D two-DOF structural model are as follows:where the symbols are defined in the nomenclature.

Nondimensionalization was applied to the time $\tau ={\omega }_{\theta }t$ and the heave displacement $\tilde{h}=2h/c$. By substituting the dimensional parameters, $I_{\theta }=M_{\theta }{b}^2{r}_a^2$; $D_{\theta }=2{\xi }_{\theta }{(I_{\theta }{K}_{\theta })}^{1/2}$; $D_h=2{\xi }_h{(M_h{K}_h)}^{1/2}$; $K_{\theta }=M_{\theta }{\omega }_{\theta }^2{b}^2{r}_a^2$; and $K_h=M_h{\omega }_h^2$, for the test cases with $M_h=M_{\theta }$ and $D_h=D_{\theta }=0$ as used by Alonso and Jameson (1994), Eqs. (1) and (2) are simplified aswhere ${\xi }_{\theta }$ and ${\xi }_h$ represent pitch and heave damping ratio, $D_{\theta }/2{(I_{\theta }{K}_{\theta })}^{1/2}$ and $D_h/2{(M_h{K}_h)}^{1/2}$, respectively; and ${\omega }_{\theta }$ and ${\omega }_h$ are the uncoupled structural natural frequencies in pitch and heave, respectively ($\mathrm{rad}/\mathrm{s}$).

In the discretization of Eqs. (3) and (4), the first derivative was calculated using a three-point backward differencing method, and the second derivative was calculated using a five-point backward stencil. Details are provided in Yuan et al. (2016).

The 3D CSD model comprises an FEM model for the solution of linear elasticity equations, and was constructed using Elmer’s FEM software package (Elmer 2019). The dynamic equation for elastic deformation of solids is written aswhere $\rho$ = density of the structure; $\overrightarrow{\mathbf{x}}$ = displacement field; $\overrightarrow{\mathbf{f}}$ = given volume force; and $\mathbf{\tau }$ = stress tensor.

In this study, a flexible 3D AGARD 445.6 wing was selected as the test case. It was modeled as a rectangular cross-sectional cantilever beam with planform dimensions (chord length, taper, swept) the same as that of the wing. Note that even though the AGARD 445.6 wing was homogenous and made of laminated mahogany, a rectangular cross-sectional wing assumption was made to ease the addition of future underwing stores and the subsequent store sensitivity analysis. Future work will be aimed at improving the structural model of the wing and the attachment of underwing stores.

The unstructured mesh of the wing geometry is comprised of second-order tetrahedron elements generated using GMSH software (Geuzaine and Remacle 2009). Since the ultimate goal of this work was to apply the FSI coupling and the subsequent sensitivity analyses to complex wing and pylon/store configurations, the wing was meshed using unstructured tetrahedron elements. Structured hex mesh could be more efficient when modeling simple geometries like rectangular beam. The second-order tetrahedron elements correspond to 3-node line elements in the tetrahedron wherein the field variable is interpolated between end nodes using a quadratic shape function. This type of adjustment is referred to as p-refinement in FEM, resulting in better conditioning and faster convergence of the FEM solution.

In this study, the linear elasticity solver in Elmer was used. It computes the structural displacement field from the Navier equations that correspond to linear theory of elastic deformation of solids. The Navier equation [Eq. (5)] was discretized in a fully implicit manner with second-order accuracy in time. The linear solver was validated by comparing eigenvalues and mode shapes for the weakened Model 3 of the AGARD 445.6 wing reported by Yates (1987). A convergence study was performed to investigate the effect of mesh density on the eigenvalues and the time history of the wing response. The mesh density chosen for the CFD-CSD coupling provided an insignificant change in eigenvalues with an increase in the mesh density.

In coupled CFD-CSD simulations using the partitioned approaches for FSI systems, the aerodynamic and structural meshes are normally separated. Proper data transfer mechanism is required to transfer the loads from the aerodynamic grid to the structural grid and the displacement computed in the structural grid to the aerodynamic grid. Similar to the aeroelastic analysis software ZAERO as described in its manual (ZONA Technology 2011), spline matrices were used in this study to relate and interpolate the displacements computed at the structural finite element grid to aerodynamic points. The spline matrix takes the following form:where $\mathbf{G}$ is the total spline matrix relating $\overrightarrow{\mathbf{h}}$ to $\overrightarrow{\mathbf{x}}$, with $\overrightarrow{\mathbf{h}}$ as the interpolated displacement vector at the aerodynamic points and $\overrightarrow{\mathbf{x}}$ as the displacement vector defined at the structural finite element grid points. Based on the principle of virtual work, the force transfer from the aerodynamic grid to the structural finite element grid are then determined by the transpose of the matrix $\mathbf{G}$:

An interpolation function was defined for the general spline problem, which is composed of a linear polynomial and a radial basis

The RBF $K_i(x,y,z)$ defines the radial interactions with the source points and the polynomial corrector $a_0+a_{1}x+a_{2}y+a_{3}z$ guarantees compatibility for rigid modes. An RBF is a real-valued function whose value depends only on the distance from the origin ($r$), so that $K_i(x,y,z)=K(r_i)$, where $r_i^2={(x-x_i)}^2+{(y-y_i)}^2+{(z-z_i)}^2$. An RBF exists if the coefficient $F_i$ and the weight of the polynomial can be found such that the polynomial terms gives zero contributions at the source points and the desired function values are obtained at source points (Cella and Biancolini 2012). This gives

Zero sum of $F_i$ in Eq. (9) can be recognized as the discrete force equilibrium equation and the others are discrete moment equilibrium equations (ZONA Technology 2011). Combining Eqs. (8) and (9) gives the following matrix system:

Built on the known values at the source points, the matrix system can be defined and the inverse matrix $[\mathbf{C}{]}^{-1}$ can be used to calculate the coefficients $F_i$ of the RBF and the coefficients $a_{0–3}$ of the linear polynomial in Eq. (8), resulting in the interpolation matrix $\mathbf{G}$ in Eq. (6). The matrix $\mathbf{G}$ was pregenerated and stored at the beginning in the simulation. A number of RBFs were implemented in the current coupling module, including polyharmonic spline $K(r)={|r|}^n$ with $n$ odd, thin plate spline $K(r)={|r|}^n\ln |r|$ with $n$ even, multiquadric $K(r)=\sqrt{1+r^2}$, inverse multiquadric $K(r)=1/\sqrt{1+r^2}$, inverse quadric $K(r)=1/(1+r^2)$, and Gaussian $K(r)=e^{-r^2}$ RBFs, as listed by Cella and Biancolini (2012). The constant “1” in the quadric splines can be varied to control the shape of the basis functions. de Boer et al. (2007) used a value of ${10}^{-6}$ in the multiquadric spline $K(r)=\sqrt{{10}^{-6}+r^2}$ in their study, which is close to the simple polyharmonic spline $K(r)=|r|$ used in this study.

A mesh deformation algorithm for CFD is required to deform the flow mesh based on the boundary displacements due to the interaction with the flexible structure. Since the aforementioned spline approach is valid for both volume mesh smoothing and surface shape changes (Cella and Biancolini 2012), the dynamic CFD mesh was updated using an RBF-morphed mesh without successive calling for a grid generator. In this study, the displacement interpolation using the RBFs was conducted for all nodes of the flow mesh. Fig. 2 demonstrates the mesh-morphing capability. The grid, originally generated using a grid generator for the SD7037 airfoil at $\alpha =3°$, was morphed to the same airfoil at $\alpha =9°$.

In Figs. 2(a and b), the shape of the SD7037 airfoil and its original mesh are depicted with thin lines in black. The shape of the 9° airfoil is shown with thick lines in red. The morphed mesh is plotted with thin lines in blue. The 9° airfoil was obtained by rotating the airfoil about its quarter chord. The original mesh was a C-type mesh with $673\times 97$ grid points. Since the 2D mesh has a small size, all 864 boundary points at farfield and all 480 points on the airfoil were used as source points to control the RBF mesh morphing. Fig. 2(c) shows the grid points on the airfoil surface, where the resulting grid points of the morphed mesh (with thin solid lines with filled circles in blue) overlapped with the source points (with thick dashed lines with open circles in red). For large-sized mesh or complex 3D wings as discussed later, one can select some of the points on the boundaries as source points. Both the selected and the rest points on the boundaries can be used to check the accuracy and correctness of the mesh morphing. In this study, the farfield points remained unchanged. The airfoil surface change, i.e., the deflection of the target airfoil shape presented by the thick dashed line with open circles in red from the original airfoil depicted by the thick solid line with open circles in black, together with the farfield points, were used as the known displacements $w$ to build the RBFs and define the interpolation coefficients ($a_{0–2}$ and $F_i$) in Eq. (8), thus generating the interpolation matrix $\mathbf{G}$ in Eq. (6). Since both farfield and wall boundaries were considered to build source points, the entire domain of the mesh was morphed or smoothed in this study. The mesh morphing was conducted using an in-house-developed mesh morphing module, which was integrated into the CFD-CSD coupling interface. Although the 2D mesh morphing process showed success in this study, it may deteriorate the mesh quality for large deflections or complex 3D deformations, in particular the mesh orthogonality as discussed later.

Depending on the choice of the RBF, the resulting coefficient matrix [C] in Eq. (10) is usually not diagonally dominant or may be ill-conditioned. When the number of source points is large, inverting this type of stiff matrix system requires significant computational effort. The resulting inverse matrix ${[\mathbf{C}]}^{-1}$ may be inaccurate causing nonphysical grid deformation as demonstrated by Yuan et al. (2016). Nevertheless, as far as products $[\mathbf{C}]{[\mathbf{C}]}^{-1}$ and ${[\mathbf{C}]}^{-1}[\mathbf{C}]$ are ensured to be an identity matrix, the source points remain unmoved during mesh morphing, which is evidenced in Fig. 2(c). Having tested the available RBFs, the simple polynomial spline RBF $K(r)=|r|$ showed the most robust behavior and delivered the best results. Consequently, the simple polyharmonic spline was used for the RBF morphing in this study. Changing the morphing domain or varying the constant (instead of 1) in the quadric functions may improve the results using other RBFs, which was not further investigated in this study.

The CFD and CSD codes were tightly coupled using a fully implicit scheme; namely, the CFD and CSD solvers exchanged data for each pseudo-time step in the simulations. Fig. 3 demonstrates the tightly coupled partitioned approach. The process was as follows:

The flow solver was first validated by simulating the flow past a forced pitching NACA 64A010 airfoil. This included the validation of the numerical discretization scheme on moving grid and the geometric conservation law implemented in the FLOWer CFD solver. The NACA 64A010 airfoil was sinusoidally oscillating in pitch about its quarter chord. The pitch angle of the airfoil was prescribed as a function of time, $\theta =0°+1.01°\cos (2\pi \mathrm{ft})$. The primary parameters used in the computation were defined according to the experiments reported by Davis (1982): chord length $c=50\mathrm{cm}$; freestream Mach number $M_{\infty }=0.796$; Reynolds number based on the airfoil chord length and freestream velocity ${Re}_c=12.56\times {10}^6$; and frequency $f=34.4\mathrm{Hz}$.

URANS simulations with Menter’s SST turbulence model were conducted for this transonic study. Simulations were performed on both an O mesh and a C mesh. The farfield of the grids was located at $15c$ from the airfoil as recommended by Wang and Zha (2010). A further increase of the farfield distance in their computations for this test case did not improve the solution. The maximum nondimensional distance ($y^{+}$) of the first grid point to the airfoil surface was $\sim 1$ in the current simulations. Wang and Zha used an O mesh with $281\times 66$ grid points. In this study, the O meshes had $721\times 145$ and $481\times 97$ grid points while the C mesh had $673\times 97$ points to ensure that the pressure gradient at the leading edge was well resolved and the shock movements at the middle of the airfoil were well captured. Simulations on the $721\times 145$ O mesh did not improve the predictions of the aerodynamic coefficients when compared with those on the $481\times 97$ O mesh. Therefore, only the $481\times 97$ O mesh and the $673\times 97$ C mesh are considered subsequently. Both had the same grid distributions on the airfoil surface with 481 grid points, but the C mesh was expected to be more suitable for capturing or modeling the shedding vortices in the wake. Since the 2D airfoil was treated as a rigid body, the mesh oscillated as a solid body without the applications of the RBF in the 2D simulations.

The time-accurate simulations were started from a stationary flow on the O mesh and from a steady solution on the C mesh. In general, the aerodynamic coefficients of the fourth cycle reproduced the ones of the third cycle. A dual time stepping approach was used in the fully implicit second-order time integration. A time step and pseudo-time dependence study was performed on the O mesh. The time step was fixed during the time-accurate simulations. Doubling the time steps from 240 to 480 for each pitching cycle did not change the aerodynamic coefficients; therefore, the results obtained using 240 time steps for each pitching cycle are shown in this paper. Normally, several hundred pseudo-time steps were required to reduce the convergence residuals by 3–4 decades for each time step. In this study, the discrepancies of the aerodynamic coefficients obtained using 160 and 640 pseudo-time steps were marginal; therefore, 160 pseudo-time steps were then used for each real time step, which resulted in a reduction in the residuals of about 1.5 orders of magnitude. This is reasonable, because, in most time instances, the flow field from the preceding time step already provides a good initial condition for the new time step, and therefore the initial residual is already low. As mentioned by Liu et al. (2001), further reduction of the residuals by 3–4 orders of magnitude was neither practical nor necessary.

Figs. 4 and 5 compare the computed results and original experimental data. For comparison, the computed results obtained by Wang and Zha (2010) using a fifth-order weighted, essentially nonoscillatory (WENO) scheme are also plotted in Fig. 4. When compared with the WENO results, the C mesh showed slightly better aerodynamic coefficient slope than the O mesh in this study. Therefore, only results computed from the current C mesh are given in Figs. 4 and 5. The plot for the pitching moment coefficient in Fig. 4 showed obvious discrepancies between numerical and experimental results. Wang and Zha believed that this large discrepancy was due to the experimental uncertainties and the inadequacy of the turbulence models, as they used the Spalart-Allmaras one-equation turbulence model in their study. Bendiksen (2011) noted the same discrepancy and attributed it to the turbulence model as well. We believe that the discrepancies were probably caused by the uncertainties (as discussed later) in both the experimental and numerical simulations due to aerodynamic nonlinearities caused by shock-boundary layer interaction. Our experience with low-Reynolds-number studies (e.g., Yuan et al. 2013) and with high angle-of-attack unsteady aerodynamics showed that aerodynamic nonlinearities manifest themselves more evidently in the pitching moment than in the lift coefficient. As stated by Bendiksen (2011), this remains an unsolved problem and needs further investigation.

Fig. 5 depicts the distributions of the instantaneous pressure coefficients as a function of pitch angles when the airfoil was pitching down. The computed pressure coefficients upstream of the shocks were underpredicted when compared against the experimental data, with a worst case on the lower surface at $x/c=0.2$, where $Cp\approx -0.38$ from CFD versus $-0.28$ from the experiment. However, the shock wave position and the pressure distribution after the shock were well predicted. No flow separation was observed in the present study.

It should be noted that the steady pressure coefficient of the reduced experimental data used in the aforementioned comparisons was corrected at one point on the airfoil lower surface after a critical scrutiny in this study. The measured instantaneous pressures on the upper surface of the oscillatory pitching case were recorded in the experimental report (Davis 1982) whereas it was not the case for those of the lower surface, and only reduced frequency pressure data are available for the lower surface. The steady pressure is needed to reproduce those instantaneous pressures by combining the reduced frequency pressure data with the fundamental (steady) ones reported by Davis (1982). Considering the symmetric configuration of the airfoil at $\alpha =0°$, the value of steady pressure coefficient $Cp=-0.777$ at $x/c=0.49$ on the lower surface was replaced by $Cp=-0.594$, the value of its counterpart on the upper surface of the experiments. Plots a and c in Fig. 6 demonstrate the correction, where the correction brought pressure distributions overlapping on the lower and upper surfaces. This correction was supported by Chen and Zha’s (2005) computational work as shown in Fig. 6(b). As a reference, Fig. 7 shows the uncorrected pressure distribution of the oscillating airfoil at $t/T=0.5$ and $\theta =-1.01°$, corresponding to its corrected value shown in Fig. 5(e). When compared with its counterpart at $t/T=0.5$ and $\theta =1.01°$ as shown in Fig. 5(a), the uncorrected pressure distributions on the upper and lower surfaces obtained from the reduced frequency pressure data shown in Fig. 7 were not fully symmetric by phase in time as shown in Fig. 5(e) for the symmetric airfoil in the symmetric oscillations about $\alpha =0°$. It is expected that the discrepancy in the aerodynamic pitching moments between the experimental and computed results shown in Fig. 4 would be reduced as the pitching moment of the experiments was calculated by integrating the pressures on the airfoil surfaces. The off values of the aforementioned pressure data were probably caused by tiny imperfections of lower surface or flow angularity of the wind tunnel, confirming high sensitivity in this transonic flow regime and difficulties for both experimental and numerical simulations. This correction practice also shows the complementary value of CFD for interpretation or uncertainty quantification of the experimental data. This is not the first time that CFD helped interpret or correct experimental data. Quadrio and Luchini (2002) and Chung et al. (2002) reported in their computational work on an annular pipe flow that integration of the velocity profiles experimentally measured in the radial direction did not yield a value of unity as documented by Nouri et al. (1994); therefore, the experimental data were rescaled with the bulk velocity.

Coupled CFD-CSD simulations were first performed for a rigid airfoil in two-DOF pitch and heave motion. The two-DOF airfoil model was first introduced by Isogai (1979) based on a 3D swept-back wing with a NACA 64A010 cross section. A number of researchers have investigated this model numerically, including Alonso and Jameson (1994), Liu et al. (2001), Isogai (1981), Bohbot and Darracq (2001), Prananta et al. (1998), Chen and Zha (2005), and Kousen and Bendiksen (1988). These studies and their numerical results are used as comparison references for the present study. Fig. 8 demonstrates a schematic of the binary system of the airfoil. According to Isogai, the structural parameters used in this model are the location of the elastic axis (normalized by semichord, positive aft midchord) $a=-2.0$, static imbalance $x_{\theta }=1.8$, ratio of the uncoupled structural natural frequencies ${\omega }_h/{\omega }_{\theta }=1$, radius of gyration (normalized by semichord) $r_a^2=3.48$, and airfoil mass ratio $\mu =60$. Note that the value of $a=-2.0$ refers to an elastic axis located upstream of the airfoil leading edge, as to mimic the dynamics of an outboard section of a swept-back wing (Alonso and Jameson 1994). The wind-off coupled frequencies are $0.713{\omega }_{\theta }$ and $5.34{\omega }_{\theta }$ for Mode 1 and Mode 2, respectively (Isogai 1979).

Both Euler and URANS simulations were carried out for the flow past the airfoil. In the URANS simulations, the Reynolds number was set to $1.256\times {10}^7$. To accelerate the computations, the coupled CFD-CSD simulations were started from the computed results of the airfoil in a prescribed pitching motion. The forced pitching motion was scheduled based on the natural structural frequency and its simulation was performed for three cycles at a certain Mach number (e.g., $M=0.825$). After the coupled CFD-CSD simulations for a series of air speeds, the critical points on the transonic flutter boundary at a given Mach number were then searched. The flutter boundary was calculated by determining the critical points where the neutral response occurs. In order to classify damped, neutral, and divergent responses of the airfoil, the speed index was used when the Mach number was fixed. To match the local value of the speed index, the freestream temperature was adjusted in the simulations, which was equivalent to adjusting the freestream pressure.

Fig. 9 shows the computed time history of pitch and heave response for the Isogai wing model for $M_{\infty }=0.825$, corresponding to the damped, neutral, and divergent oscillatory responses, respectively. When the speed index was below the critical value on the flutter boundary, both pitch and heave displacements decayed, indicating a damped response as shown in Fig. 9(a). In Fig. 9(b), the neutral response appeared periodically with a well-defined frequency because the speed index was close to the critical value. Once the speed index was beyond the critical value, a divergent oscillation occurred, as shown in Fig. 9(c). The nondimensional flutter frequency was calculated to be ${\omega }_f/{\omega }_{\theta }=0.85$ in this study. It lies well between 0.78 (Isogai 1981) and 0.88 (Prananta et al. 1998), where the numbers are digitalized readings from Fig. 12 in Prananta et al. (1998). In comparison to the wind-off modal frequencies, it suggests that the flutter mode was Mode 1, as discussed by Prananta et al. or by Weatherill and Ehlers (1983). However, some ambiguity remains concerning the nature of the flutter mechanism. Starting with Isogai (1979), most researchers describe Mode 1 as heave (or plunge, or heave dominated, or bending), whereas Weatherill and Ehlers (1983) clearly showed that a single DOF simulation in pitch led to a very similar flutter mechanism both in flutter speed index and frequency.

Table 1 lists the computed critical speed indices from the Euler and URANS simulations. The computed speed indices compare reasonably with the results from other researchers. The viscous results of the critical speed index computed on both O and C meshes (0.685 and 0.72, respectively) lay well among the URANS results (0.615–0.82) from other researchers. The inviscid results of the critical speed index computed from the O mesh (0.57) was slightly underpredicted and the one calculated on the C mesh (0.605) showed a slightly better agreement when compared against the solutions from other researchers (0.60–0.64). This is reasonable as the grid lines on the C mesh have a better alignment and orthogonality.

Numerical convergence was monitored to ensure that the coupled CFD-CSD simulations were performed in a time-accurate manner. Fig. 10 demonstrates the convergence history of a time step after the flow reached a state of neutral oscillatory responses. As discussed in the previous section, 160 pseudo-time steps were limited for each real time step. Fig. 10 confirms a reduction in residuals of about two orders of magnitude within 160 iterations (pseudo-time steps), which was the same as discussed for the prescribed pitching case. In the figure, the residual of the continuity equation was normalized by that of the start flow field, which was initialized using freestream quantities. A closer look showed that after eight pseudo-time steps with a residual drop by 0.5–1 order of magnitude, changes in aerodynamic coefficients were marginal in this study. Therefore, 10–20 pseudo-time steps are necessary for coupled CFD-CSD simulations used for flutter analysis.

The computations were repeated for several Mach numbers and a flutter boundary was determined as shown in Fig. 11. Several runs at different speed indices were required to identify each point on the boundary. The identification of the boundary was relatively straightforward for the lower branch, displayed as symbols with a connecting line. This is exemplified in Fig. 9, where the value of the speed index that yielded a zero damping response is clear. On the other hand, the boundary defined by the nonconnected symbols is not trivial. In these cases, the dynamic behavior is indicative of strong nonlinear effects, whereby the boundary is no longer defined by the loss of stability of the equilibrium point (as seen in linear flutter) but as the loss of stability of limit cycle oscillation.

Fig. 11 compares the computed flutter boundaries from the current study with previous references (Isogai 1981; Bohbot and Darracq 2001; Prananta et al. 1998). Known as the transonic dip due to compressibility effects, there was a drop in the flutter boundary of the aeroelastic system placed in a transonic flow. As can be seen, the transonic dip phenomenon was well predicted in the present CFD-CSD simulations. The bottom of the dip was located at a speed index of $V_f=0.54$ in the current Euler simulation. As shown in the figure, the transonic small disturbance (TSD) approach predicted a higher speed index in the low-Mach range ($M_{\infty }\sim 0.75$), possibly due to the inadequate prediction of the phase lag in the shock wave motion (Alonso and Jameson 1994). The TSD and the Euler solutions seemed to agree more closely at $M_{\infty }\sim 1$, probably due to the fact that the TSD is an exact approximation in the limit $M_{\infty }\to 1$. Disagreement of the present inviscid solution with Prananta et al.’s Euler results was attributed to the different grid densities. Prananta et al. (1998) used a mesh with $140\times 60$ grid points in their simulations.

As in the studies reported by Prananta et al. (1998) and Alonso and Jameson (1994), multiple flutter points were observed for a range of Mach numbers in the neighborhood of $M_{\infty }=0.85–0.875$ for the inviscid simulations, which is evidenced of a strongly nonlinear problem. For higher Mach numbers, flutter points existed but occurred at higher speed indices. In fact, the flutter mode was then the second mode of vibration with a much higher frequency as evidenced in Fig. 12, where ${\omega }_f=5.03{\omega }_{\theta }$. In Fig. 12, the initial low-frequency oscillations eventually converged to high-frequency oscillations. When compared with Fig. 9, in addition to the high frequency, the phase shift is obvious between the heave and pitch oscillations. The convergence to the higher frequency oscillations with phase shift will be further discussed later.

As mentioned previously, LCOs were observed in the Euler simulations as reported in the public literature, e.g., Alonso and Jameson (1994), Liu et al. (2001), Isogai (1979), and Kousen and Bendiksen (1988). In this study, as seen in Fig. 13, the vibrations initially diverged but eventually reached steady-state limit-cycle oscillations due to nonlinearities in the aerodynamics after the long time computations. This exercise confirmed the ability of the current CFD-CSD coupling algorithm. The nondimensional LCO frequency was calculated to be ${\omega }_f/{\omega }_{\theta }=1.22$ for $M_{\infty }=0.825$ in this study, comparable to the digitalized reading of ${\omega }_f/{\omega }_{\theta }\approx 1.34$ for $M_{\infty }=0.75$ from Liu et al. (2001). No phase shift was observed between the heave and pitch oscillations in the LCO. Alonso and Jameson (1994) had computed the $M_{\infty }=0.75$ case as Liu et al. (2001) but showed a very low frequency of ${\omega }_f/{\omega }_{\theta }\approx 0.15$ with a phase shift between the heave and pitch. It is believed that there was an improper normalization in their postprocessing, and their phase shift could be attributed to the inconsistent sign definitions for the pitch and/or heave motion as the phase shift was close to 180°.

The appearance of a second stable LCO branch at a higher frequency for higher Mach numbers is a consistent observation found in the open literature. In our simulations, the solution was started from a prescribed pitching motion for at least three oscillation cycles at a frequency equal to the uncoupled structural frequency in pitch, which corresponds to the first LCO branch (Fig. 13). The aeroelastic system was then freed to respond freely. At a Mach number of 0.9, only the high-frequency LCO exists. Hence, the time response evolved from the low-frequency transients to a high-frequency steady-state oscillation regime (Fig. 12).

Computations of the viscous flutter boundary were performed using Menter’s SST model. The predicted flutter boundary is also depicted in Fig. 11. As can be seen in the figure, the bottom of the dip was located at a speed index of $V_f=0.72$ ($M_{\infty }=0.825$) from the present URANS computations, while it was $V_f=0.54$ ($M_{\infty }=0.85$) in the Euler simulations. This indicates that viscosity reduced the Mach number for the transonic dip. Although the computed viscous flutter boundary compared well with the results from other researchers, the discrepancies among the solutions were visible, especially at high speed indices. It is believed that the discrepancies were mainly due to the use of different numerical schemes because the solution using the Spalart-Allmaras one-equation turbulence model did not show a big difference from the SST results in this study. Furthermore, the viscous and inviscid simulations are in relatively good agreement up until the bottom of the transonic dip. For higher Mach numbers, it is clear that viscosity plays an important role, and the inviscid calculations are no longer reliable. In particular, the S shape of the instability boundary in the transition dip region predicted by the Euler simulations no longer exist when viscosity is introduced.

The cantilevered AGARD 445.6 wing was chosen to validate the coupled CFD-CSD solver and investigate the 3D aeroelastic phenomena. The aforementioned 2D simulations from the previous sections were used as reference for the setup of 3D computations reported in this section. The AGARD 445.6 wing model has a 2.5-ft semispan and consists of a symmetric NACA 65A004 airfoil, with a sweep angle of 45° and a taper ratio of 0.66 (Yates 1987). The chord lengths are 1.833 and 1.208 ft at the wing root and tip, respectively. In this study, the leading edge of the root was set at the origin of the coordinate system. The computational domain was discretized using multiblock structured grids. Based on the previous 2D simulations, the 3D computational domain around the wing consisted of $673\times 97\times 33$ grid points. To resolve the wing tip flow, the flow domain was extended by one semispan. The extended computational domain consisted of $673\times 97\times 33$ plus $241\times 17\times 33$ grid points, with a farfield located at 15 local airfoil chords from the wing surface. The computational domain was partitioned into 29 blocks using parallel computing to expedite the coupled simulations. Fig. 14 demonstrates the deflected 3D AGARD 445.6 wing with its sectional meshes at the root and tip. The 2D meshes are sectional cuts from the 3D mesh after morphing.

An RBF morphing technique was used for the 3D mesh deformation. The farfield boundaries were set as static while the wing surface was set as a flexible deflection boundary to be defined by the 3D structural model as discussed subsequently. It was observed in a previous work on prescribed-motion computations of the AGARD wing “weakened Model 3” that the difficulties in mesh deformation occurred mainly due to the second torsion mode (Yuan et al. 2016). To demonstrate the severity of the mesh deformation, the 3D mesh shown in Fig. 14 was morphed based on the second torsional mode (Mode 4) according to Yates (1987). The deflection amplitude was defined such that the maximum deflection, which is at the trailing edge of the tip for the fourth mode, was as large as 20% of the length of the semispan (i.e., 40% of the chord length at the wing tip). As the root was fixed mimicking in the experiment, the grid orthogonality at the root cross section was nearly unchanged when compared to the original base mesh [Fig. 14(c)]. However, the grid orthogonality was severely deteriorated at the wing tip [Fig. 14(d)]. The angle between the cell face normal and the line connecting the cell centers may become close to or even larger than 45°, which could affect numerical accuracy. In this situation, allowing the farfield boundary to adapt to some extent according to the wing body deflection would improve numerical accuracy. Because the deflections in the current 3D CFD-CSD computations were generally gentler than the case demonstrated in Fig. 14, which was true at least for the neutral and damped oscillations, the farfield boundary adaption was not considered in this study. Note that the resolution in the boundary layer of the morphed mesh remained comparable to the original base mesh as the morphing process is actually an interpolation depending on distances between points, meaning short distance remains short.

The “weakened Model 3” of the AGARD wing is a well-known test case considered as an AGARD standard aeroelastic configuration for flutter calculations (Yates 1987). In this study, the AGARD 445.6 wing was modeled structurally as a 2D plate or a 3D elastic solid, as demonstrated in Fig. 15. The 2D plate model is based on the shear deformable model of Reissner (1945), Mindlin (1951), and Reddy (2006). The 3D FEM model was based on the linear theory of the elastic deformations of solids, namely, the Navier equations, wherein the finite element mesh is discretized using second-order tetrahedron elements. Since the AGARD wing was made of mahogany, the wing material was considered to be orthotropic wherein the primary material axis was aligned along the midchord line. The Young’s Modulus along various axes and Poisson’s ratio were slightly adjusted from the reported values to align the first four eigenvalues in close agreement with the published experimental and numerical values. Also, the density of the wing was adjusted to make sure the mass of the wing is the same as in the experiments. The cross-sectional thickness can be constant or varying along the span as a function of the sectional chord length. For the 3D solid model with constant thickness demonstrated in this study, the wing boundary was discretized with 1,716 grid points.

The modal characteristics of the AGARD 445.6 wing were validated using both time-marching and eigenanalyses. Fig. 16 shows the time-marching solution for the trailing edge at the wing tip without any aerodynamic considerations. Table 2 compares the first four modal frequencies obtained using eigenvalue analysis on the CSD model with the analytical results from Yates (1987). The eigenvalue analysis was performed using linear elasticity equations to model the cantilevered wing beam, wherein Elmer software was used to perform direct solve on the eigenequation obtained following finite element discretization. Fig. 17 compares the mode shapes of the weakened Model 3 of the AGARD 445.6 wing calculated in the eigenanalysis for the 3D model (constant thickness) with those by Yates (1987). As can be seen in the table and figures, the frequencies and mode shapes of the AGARD wing test bed are in reasonably good agreement with the published results.

The advection upstream splitting method was used in CFD because it resulted in slightly better results when compared with the central differencing scheme. As described previously, the computational mesh had about 4.4 million grid points, which is much denser than the grid with 176,601 points used in Liu et al.’s work (2001). Fig. 18 demonstrates the AGARD 445.6 wing aerodynamic surface interacting with the structural cantilevered wing beam model. In CSD, the 3D wing was modeled as a rectangular cross-sectional cantilever beam with planform dimensions (chord length, taper, swept) the same as that of the wing. As a result, the shown structural model was not the same as a traditional wing box constructed within the real wing surface. This was chosen to more intentionally prove the concept that the structural model can be more “arbitrarily” designed as far as it closely follows the aerodynamic configuration. It was driven by the belief that the rectangular cross-sectional assumption would ease the addition of future underwing stores and the subsequent store sensitivity analysis. In fact, as far as the data transfer between CFD and CSD conserves the force, moment, and energy, this configuration works well as discussed subsequently. Nevertheless, it is straight-forward to shorten the sectional length to align the structural model within the aerodynamic configuration if wanted.

As mentioned, for the mesh deformation using the RBF morphing technique, the farfield boundaries were set as static while the wing surface was set as a flexible deflection boundary. The coarsening ratio or reduction factor to the control (source) points for the static boundary was set to 16. This means that one of every four points in the two surface coordinates were selected as the source points, resulting in 1,944 source points on the farfield boundaries. The selection of the source points on the wing surface plays an important role in applying RBF to the mesh morphing. The number of source points decides the size of the transformation matrix $\mathbf{G}$ as defined in Eq. (6), which is an inversed matrix of RBF $K_i(x,y,z)$. If the number of source points is too large, the inversion of the matrix would limit the effectiveness of the RBF mesh morphing. For very complex configurations, one may need to use multiple blocks or multiple interfaces to build the transform matrices. As a reference, $\sim \mathrm{4,000}$ source points on the wing surface were selected to test the effectiveness of the RBF morphing strategy in this study, with a coarsening ratio of four for the streamwise and one for the spanwise directions or two for both streamwise and spanwise directions, respectively, resulting in $\sim \mathrm{6,000}$ source (control) points in total for the interpolation matrix. The inversion of the matrix required $\sim 30\min$ of the CPU time on an Intel Xeon E5-2650 (2.6Hz) computer. Note that this matrix inversion is a one-time job conducted at the beginning of the coupled CFD-CSD simulations. It can also be pregenerated and saved for other similar computations. Since the source points were selected directly on the wing surface, the deflection must be first transferred from the structural surface to the aerodynamic surface before the RBF mesh morphing. Another effective way is using the structural displacements directly as source points for the RBF mesh morphing (van Zuijlen et al. 2007). This is the approach that was mainly applied in this study. As mentioned, the structural surface has 1,716 grid points, resulting in a transformation matrix with 3,660 source points in total. A CPU time of $<0.1\mathrm{s}$ was used for the CSD part and its data transfer with CFD, which was negligible. However, $\sim 10\mathrm{s}$ of CPU time were needed for deforming the whole computational domain. For more complex configurations, a zonal approach for the RBF mesh morphing and combined explicit-implicit CFD-CSD coupling mechanism may be considered or developed.

Smith (2012) pointed out that both the balance of the forces and the equivalence of the total work between the CFD and CSD methods are necessary. As discussed previously, Eq. (9) serves as the discrete force and moment equilibrium equations while the work conservation is ensured by Eq. (7). To verify the discretization of the equations, the integrated aerodynamic forces, moments, and real work on the interface were monitored during the CFD-CSD simulations, and the results are exemplified in Table 3. Table 3 compares the aerodynamic forces and moments on the AGARD 445.6 wing surface and those on the CSD model wing surface at representative time instances for two test cases for $M_{\infty }=0.901$. The aerodynamic force on the wing surface was calculated from the CFD solver by integrating the pressure distributions while the shear part was considered marginal. The pressure forces obtained in the CFD simulations were then transformed to the surface of the CSD model using the spline matrix defined in Eq. (7), which is then used to compute the structural FEM solution. The respective forces were used to calculate the aerodynamic moments about the reference axes. The aerodynamic forces and moments were then summed in three coordinates. The real work was calculated by the respective forces along the deflections. Table 3 confirms the force, moment, and work (energy) conservations in the aerodynamic load transferal on the CFD-CSD interface, indicating effectiveness of the aforementioned interface coupling scheme applied to the AGARD 445.6 wing aerodynamic surface interacting with a structural cantilevered wing model.

The CFD-CSD simulations were performed for the AGARD 445.6 wing at 0 degrees of angle of attack using a number of structural models. For the sake of brevity, this paper is focused on the 3D model with a constant thickness along the span, and only the vertical component of the pressure force was considered. Detailed comparisons among various structural models including the effects of the underwing stores are out of the scope of the present paper. In analyzing the (flutter) stability boundary, it was observed that the spanwise and chordwise components of the aerodynamic loads had a minimal effect on the temporal response of the wing due to small deformations, and therefore were ignored in this paper. Fig. 19 depicts the instantaneous pressure distribution at the midspan and near the tip obtained from one of the simulations, indicating the occurrence of the shocks in the transonic regime. Fig. 20 illustrates the computed deflections at the wing tip and the lift coefficients of the elastic AGARD 445.6 wing at $M_{\infty }=0.901$. Results shown in Fig. 20 were obtained from the fully coupled 3D CFD-CSD simulations. Simulations were started from a static condition, followed by a motion schedule prescribed according to Mode 1, or combined with Mode 2, as listed in Yates (1987). Modes 1 and 2 were chosen to accelerate the simulation convergence since Yates (1987) indicated that the flutter frequency lies between structural Modes 1 and 2. The flow field was believed to be relatively well developed after one to two cycles of the prescribed motion. Subsequently, the prescribed schedule was removed and the system became an elastic one.

The plots in Fig. 20 show the damped, neutral, and divergent responses, respectively. The flutter boundary was then computed by determining the neutrally stable or critical point when the neutral response occurred. When the speed index was below the critical value on the flutter boundary, both displacements and lift coefficient decayed, indicating a damped response as shown in Fig. 20(a). In Fig. 20(b), the neutral response appeared periodically with a well-defined frequency as the speed index was close to the critical value. Once the speed index was beyond the critical value, a divergent oscillation occurred, as shown in Fig. 20(c).

After repeating the process for several Mach numbers, a flutter boundary for the AGARD 445.6 wing was determined. Several time-marching runs were needed to compute the flutter speed index for each Mach number by manually changing the freestream temperature, thus the flight speed, until a zero damping response was found. In the current study, the upwind MAPS method (Rossow 2000) was mainly used in the transonic flow simulations. However, numerical instabilities were encountered during the simulations for the conditions with $M_{\infty }>1$. In these cases, the advection upstream splitting method (AUSM) (Kroll and Radespiel 1993), similar in principle to the van Leer flux vector splitting scheme (van Leer 1979), was applied instead. Fig. 21 compares the calculated flutter boundaries from the coupled CFD-CSD simulations with the experimental and analytical data from Yates (1987). For comparison, the computational results from Chen et al. (2007) and Liu et al. (2001) are also included. Both the trend of the flutter boundary and the transonic dip were well predicted in the present simulations when compared to the experimental data. The bottom of the dip was obtained at a speed index of $V_f=0.29$ in the current study, while it was $V_f\sim 0.30$ in the experimental measurements. This discrepancy was about 3%. The simulations did not consider structural damping, which would raise the flutter speed slightly. When compared with other numerical results, the present simulation showed a slight advantage over the other listed numerical methods for Mach numbers at and post the bottom of the transonic dip. At $M_{\infty }=1.14$, the other two CFD-CSD simulations overpredicted the flutter index while the current simulation quite accurately predicted it. It is believed that the discrepancy of the Euler results obtained by Liu et al. (2001) was mainly attributed to neglecting viscous effects in the inviscid simulations. Such discrepancies were also observed in the inviscid simulations conducted by other researchers, e.g., Silva et al. (2014) and Lee-Rausch and Batina (1995), where the flutter dip at $M_{\infty }<1$ and the rise in flutter speed index at $M_{\infty }>1$ were largely overpredicted.

As a side note, Silva et al. (2014) referred the flutter boundary approaching to $M_{\infty }\sim 1$ to compressibility flutter dip rather than the traditional transonic flutter dip, because the AGARD 445.6 wing has a thin airfoil, and they believed that the wing does not reach transonic condition until about $M_{\infty }=0.98$. In the current work, we still use the conventional term “transonic flutter dip” as we observed supersonic flows (local $M>1$) at some locations in the simulations for $M_{\infty }\ge 0.9$, similar to illustrations in Fig. 19, which needs to be further investigated.