2.1 Linearised frequency-domain solver

The LFD approach is first introduced. For a finite-volume method, the semi-discrete RANS equations can be written as

(1) $$\begin{eqnarray} \frac{{\rm d}{M(x){{w}}}}{{\rm d}t} + {R}({w}, {x}, \dot{ {x}})&= 0, \end{eqnarray}$$

with the diagonal matrix M storing the cell volumes and the residual function R depending on the vectors of fluid unknowns w, grid-point locations x and grid-point velocities $\dot{{{x}}}$ . Assuming small perturbations $({{w}}_1,{{x}}_1,\dot{{{x}}}_1)$ from a steady state (w ₀, x ₀), the variables can conveniently be separated as

$$\begin{eqnarray*} {{w}}(t) = {{w}}_0 + {{w}}_1(t), \qquad {{x}}(t) = {{x}}_0 + {{x}}_1(t), \qquad \dot{{{x}}}(t) = \dot{{{x}}}_1(t)\nonumber \end{eqnarray*}$$

and Equation (1) can be linearised around this steady state

(2) $$\begin{eqnarray} M \frac{{\rm d}{{{w}}_1}}{{\rm d}t} + {{R}}\left({{w}}_0, {{x}}_0\right) + \frac{{{{\partial R}}}}{{{{\partial w}}}} {{w}}_1 + \frac{{{{\partial R}}}}{{{{\partial x}}}} {{x}}_1 + \frac{{{{\partial R}}}}{\partial \dot{{{x}}}} \dot{{{x}}}_1 + {{w}}_0 \frac{\partial M}{{{\partial x}}} \dot{{{x}}}_1 = 0 \end{eqnarray}$$

The residual at steady state R(w ₀, x ₀) is assumed to be negligible small and Equation (2) is then transferred into the frequency domain

(3) $$\begin{eqnarray} \left[\frac{{{{\partial R}}}}{{{{\partial w}}}} + j\omega ^\ast M\right]\widehat{{{w}}} = -\left[ \frac{{{{\partial R}}}}{{{{\partial x}}}} +j\omega ^\ast \left(\frac{{{{\partial R}}}}{\partial \dot{{{x}}}} + {{W}}_0 \frac{\partial M}{{{\partial x}}}\right)\right]{{\hat{x}}} , \end{eqnarray}$$

with j denoting the imaginary unit. The reduced frequency ω* has been normalised with the chord length and the free-stream velocity

(4) $$\begin{eqnarray} \omega ^{\ast} = \frac{2\pi f\,l_{\text{ref}}}{U_{\infty}} \end{eqnarray}$$

Equation (3) relates the Fourier coefficients of a harmonic excitation ${{\hat{x}}}$ to the Fourier coefficients of the fluid unknowns $\widehat{{{w}}}$ , constituting a large, but sparse, linear system of equations. In the frequency-domain approach, the flux Jacobian matrix is evaluated analytically, while the linearisation with respect to grid motion and velocity uses central finite differences. Including the linearisation of the turbulence model is crucial to accurately capture viscous effects like shock boundary-layer interactions. Further details concerning the LFD method in the DLR-TAU code and effects of simplifications within the Jacobian $\frac{{{{R}}}}{{{{w}}}}$ can be found in Ref. 13. The accuracy of this method depends on two assumptions: a dynamically linear response and a negligible steady-state residual. Especially the latter might be violated for complex geometries at challenging flow conditions. In a preliminary study, the mean solution of a time-accurate simulation was used in as linearisation point, when the steady solver did not converge, yielding still reasonable results.^{(Reference Thormann, Nitsche and Widhalm11)}

The central part when applying an LFD method is solving the system of linear equations, corresponding to the exact order of the underlying spatial scheme, efficiently in terms of computational time and memory requirements. A generalised conjugate residual solver with deflated restarting is used, which recycles an invariant Krylov sub-space between restarts of the underlying generalised minimal residual solver.^{(Reference Xu, Timme and Badcock17)} The linear system is pre-conditioned using an incomplete lower-upper factorisation of a blended flux Jacobian matrix resulting from first- and second-order spatial discretisations.^{(Reference McCracken, Ronch, Timme and Badcock18)} Applying these methods significantly improved the robustness, providing convergence even at severe flow conditions.

2.2 Pulse excitation

A common approach to identify frequency response functions with a time-domain solver is to use a sinusoidal excitation while integrating in time until the response becomes periodic. This process is repeated for all frequencies of interest. However, if small disturbances are assumed, the aerodynamic system responds dynamically linear. The superposition principle can then be applied leading to a more efficient approach. A non-periodic time-domain function (e.g. a pulse, chirp or step) can be used to excite a broad frequency range. Hence, a frequency response function H per mode can be obtained from one single time-domain simulation when dividing the Fourier transform of the output signal by the Fourier transform of the excitation signal, e.g. for an arbitrary response ζ due to arbitrary input q,

(5) $$\begin{eqnarray} H_{\zeta ,q}(\omega ^\ast ) = \frac{\mathcal {F}(\zeta (t)-\zeta _0)}{\mathcal {F}(q(t))}, \end{eqnarray}$$

with $\mathcal {F}$ denoting the Fourier operator. Since this approach assumes a linearly responding system, it belongs to the group of time-linearised methods, while avoiding an explicit linearisation of the underlying governing equations.

While the particular shape of the excitation function is not important, three criteria should be satisfied nevertheless. First, its Fourier transform should not exhibits roots in the magnitude at frequencies within the range of interest. An example is presented in Fig. 1(a) showing the Fourier transform of four excitation functions. While the chirp and the step functions result in a nearly constant magnitude over the relevant frequency range, the magnitude of the 1-cos function shows two roots and a significant decrease with increasing frequencies. Such behaviour can be avoided if a non-symmetric polynomial is used instead. The non-symmetric polynomial used in this study takes the specific form

$$\begin{eqnarray*} q(t) = (6t^2 -15t + 10)t^3 , \quad t \in [0,1]\nonumber \end{eqnarray*}$$

for the ascending part of the pulse. The function is constructed with the boundary conditions, that the first and second derivative is zero at both endpoints of the interval. The final form is obtained by mirroring and stretching the polynomial for the descending part by a factor of three. Its time-domain representation as well as the 1-cos function are provided in Fig. 1(b). Both functions have the same compact support, while the polynomial exhibits its maximum earlier than the 1-cos function.

Figure 1.

Fourier transforms of selected excitation signals and corresponding time histories of pulse excitations.

A compact support is the second criterion. Such support is preferred to reduce the number of grid deformation and preprocessing calls while integrating in time. Preprocessing is usually required to update the dual-grid metrics, but it is relatively fast compared with the deformation. For the wing-fuselage configuration discussed below, grid deformation and pre-processing combined account for almost as much computational time per physical time step as the iterations of the non-linear flow solver, which underlines the necessity to perform as few grid deformations as possible. Thus, using a chirp function as excitation, which offers almost ideal frequency content, would lead to significant computational overheads.

In the application of CFD solvers, smooth functions are preferable. While the chirp and the two pulse functions satisfy this condition, the step function does not. In the case of a step function, the grid velocity during the first time step is

(6) $$\begin{eqnarray} \dot{x_1} \approx \frac{x_1(t+\Delta t) - x_1(t)}{\Delta t} = \frac{1}{\Delta t} \end{eqnarray}$$

and hence inverse proportional to the time-step size. When severe flow conditions demand very small computational time steps, the grid velocity can become large, causing serious convergence issues for the time integrator. Therefore, applying a step function would constrain the minimum time-step size. Consequently, the non-symmetric polynomial pulse, satisfying all criteria, is used as excitation function throughout in this paper.

Applying the pulse technique comes with the assumption that a perturbation from the steady base state is caused by the excitation only. However, the initial steady computation is usually not converged to machine precision and additional fluctuations can occur, e.g. due to reflections from the far-field. These have an effect on the computed frequency response function especially if the excitation amplitude is small such as required at flow conditions exhibiting strong shock-induced separation. The calculation of response functions is improved by performing an additional time-dependent, static simulation. The perturbation of the output is computed by subtracting both signals instead of considering the difference to the steady solution only. Although the computational cost is nearly doubled, significant improvements can be observed at small and medium frequencies for the wing-fuselage case in pre-buffet, as shown in Fig. 2. While more detail for this configuration is discussed below, improved results are observed for reduced frequencies between 1 and 4 as well as for the quasi-steady part at zero reduced frequency.

Figure 2.

Magnitude of lift frequency response for wing-fuselage configuration.

2.3 Computational fluid dynamics solver DLR-TAU

The DLR-TAU code^{(Reference Gerhold, Galle, Friedrich and Evans19,Reference Schwamborn, Gerhold and Heinrich20)} is a finite-volume Euler and Navier-Stokes solver using hybrid grids. The chosen discretisation employs the modified scheme of Jameson, Schmidt and Turkel^{(Reference Jameson, Schmidt and Turkel21)} for the mean-flow equations, while the Spalart-Allmaras^{(Reference Spalart and Allmaras22)} one-equation turbulence model is used for the eddy-viscosity closure throughout in this paper. The flow equations are marched to steady state with a lower-upper, Symmetric-Gauss-Seidel pseudo-time integration method and geometric multi-grid. Time-accurate un-steady flow solutions can be obtained following the dual time-stepping approach^{(Reference Jameson23)} combined with a second-order accurate backward differencing scheme.