A harmonic balance method for nonlinear ﬂuid structure interaction problems Computers and Structures

Over the past decade the Harmonic-Balance technique has been established as a viable alternative to direct time integration methods to predict periodic aeroelastic instabilities. This article reports the progress made in using a frequency updating procedure, based on a coupled ﬂuid-structural solver using the Harmonic-Balance formulation. In particular, this paper presents an efﬁcient implicit time-integrator that accelerates the convergence of the structural equations of motion to the ﬁnal solution. To demonstrate the proposed approached, the paper includes a detailed investigation of the impact of input parameters and exercises the method for two types of ﬂuid-structural nonlinear instabilities: transonic limit-cycle oscillations and vortex-induced vibrations.


Introduction
The ever growing capability of computing hardware and software, enabled high fidelity computational fluid dynamics (CFD) to become the primary tool for the study of fluid physics. With similar advancements in computational structural dynamics (CSD) and coupling algorithms, CFD has also been extensively applied to fluid-structure interaction problems where flow nonlinearities such as shock waves or flow separation play a dominant role. The time dependency of this type of analysis means it usually requires additional computational resources. Nevertheless, fluidstructure interactions play a safety critical role in several applications such as aircraft flutter or vortex induced vibration (VIV) frequency lock-in. Therefore, the efficient prediction of these and associated phenomena attracts much attention from the numerical methods research community.
Several alternatives to full time-domain simulations have been proposed to investigate the loss of linear stability of an aeroelastic system. Badcock et al. avoid conducting time-domain analysis by tracking the critical eigenvalue of the coupled system Jacobian matrix [1][2][3], at an equivalent cost of a few steady-state calculations; if the Jacobian of the coupled system is not available, linear reduced-order models (ROM) can be built from prescribed time domain simulations [4,5] and used to infer the system's stability at a lower cost than using the full-order model.
If nonlinearities are present, additional instabilities have been observed which pose further challenges to the development of efficient simulation tools [6]. In particular, nonlinearities such as shocks, vortex shedding or free-play can cause periodic oscillations, generally referred to as limit cycle oscillations (LCO) that are not captured by linear ROMs. Therefore, the development of nonlinear ROMs is a popular topic of research [7-10]. The system's nonlinear response can be approximated based on input-output relations using a recursive nonlinear interpolator in lieu of the full order aerodynamic model [11][12][13]; Yao and Marques applied an input-output technique employing a radial basis function (RBF) neural network together with a discrete empirical interpolation method to reconstruct the complete flow field for the prediction of LCO [14].
The aforementioned nonlinear methods rely on performing time domain simulations a priori, under specific conditions that enable a suitable nonlinear model for the problem to be built. Alternative methods have been proposed that preserve the underlying physics of the problem and allow the direct computation of the nonlinear system. Following the stability analysis described in Refs. [1,3], LCO can be predicted by model order reduction using the critical eigenbasis of the Jacobian of the coupled system. Another alternative in this class of methods is to exploit the periodicity of the problem and solve the fluid-structural system in the frequency domain, these strategies are commonly known as Time-Spectral method or high-dimensional Harmonic Balance (HB) approach. The HB method transfers the time dependent problem into a steady solution based on Fourier expansions. Originally applied in the context of CFD for turbomachinery problems [15,16], the HB method has been used in several diverse applications, including the prediction of dynamic derivatives [17,18], periodic flows found in rotorcraft [19] or wind-turbines [20]. Progress has also been made in applying the HB methodology to fluid structure interaction. For the majority of such cases, it is worth noting that the fundamental frequency of the oscillation is not known a priori and is an additional unknown. Thomas et al. developed the HB/LCO method to predict LCO based on a Newton-Raphson approach, where a HB formulation solves the fluid problem and amplitude and frequency are determined simultaneously [21]. Blanc et al. proposed a fully coupled aero-structural Harmonic Balance solver for forced motions. Ekici and Hall also developed a fully coupled fluid-structure coupled HB strategy capable of predicting LCOs for one degree-of-freedom (DOF) turbomachinery components [22]. Yao and Marques, building upon these approaches, proposed an Aeroelastic-HB (A-HB) to analyze fixed wing nonlinear aeroelastics [23] and this was further adapted to VIV lock-in by Yao and Jaiman [24].
The A-HB approach and its version developed for VIV employs an explicit scheme to resolve the structural equations of motion, together with a relaxation approach to update the fundamental frequency of the oscillation. Both these strategies limit the efficiency of the iterative scheme. This paper addresses these limitations by reformulating the A-HB method using an approximate exponential time differencing scheme and subsequent modified strategy to update the frequency. The following sections will describe the numerical formulation for the flow and structural models, this is followed by the introduction of the new coupling procedure; the final two sections of the paper present a diverse range of test cases used to critique the new A-HB method and the conclusions obtained from this work.

Fluid governing equation
The governing fluid equations used in this work are the compressible Euler or laminar Navier-Stokes equations. For timedependent problems with moving boundaries, the system of equations is solved using an arbitrary Lagrangian-Eulerian (ALE) formulation as follows: the integral form of these equations for a control volume V j with surface dS can be written as: where t is the physical time, w ¼ q; qu; qE ½ T is the vector of conserved variables, the over-bar denotes the control volume average quantities; q is the density and E is the energy, u and v g are the Cartesian flow and grid velocities vectors, respectively and n is the outward unit normal of every cell face. F i and F v correspond to the inviscid and viscous fluxes, respectively and can be written in compact form as: The inviscid flux is calculated by the AUSM þ À up flux function [25], together with the Van Albada limiter to achieve 2nd order spatial accuracy, details of this implementation can be found in Ref. [26]. The viscous stress tensor s ij and the heat flux q i are given by: where l is the coefficient for laminar viscosity, c p is the specific heat ratio for a perfect gas and Pr ¼ 0:72 is the laminar Prandtl number. The pressure is given by the ideal gas law: ð9Þ where x is the system's fundamental frequency, w hb and R hb are the fluid and residual solution at equally spaced time intervals Dt ¼ T=ð2N H þ 1Þ where T ¼ 2p=x and N H corresponds to the number of harmonics selected. The matrix D is given by: Full details of the HB implementation and derivation of matrix D can be found in Refs. [23,27]. The Eq. (9) is solved by introducing a pseudo time variable, s: This means that the time dependency of the fluid equation is converted into a steady problem. An LU-SGS [28] scheme is employed to March the Eq. (11) forward in pseudo-time.

Structural governing equations
The equations of motion for a linear structural model can be described in general terms as here M; f; K are the mass, damping and stiffness matrices of the structure, g and f are the displacement and external force applied to the structure. The latter two quantities also correspond to the output and input to the structural system. The structural equations in state space format can be described as, Similar to the fluid equations, the same HB operator D can be applied to Eq. (13), resulting in