2.1.1 Problem formulation

The implemented potential flow model consists of a three dimensional panel code based on the formulation described in Katz and Plotkin^{(Reference Katz and Plotkin29)}. In the panel method formulation, the flow is assumed to be incompressible, inviscid and irrotational. Considering these assumptions, the governing equation can be written in terms of velocity potential, Φ(x, y, z), by the Laplace’s equation:

(1) $$\begin{equation} \nabla ^2\Phi = 0. \end{equation}$$

2.1.2 Steady formulation

To solve the Laplace’s equation, two boundary conditions are prescribed:

• • impermeability condition (zero normal velocity on a body),

• • far field condition (the motion should disappear far from the body).

The Laplace’s equation (Equation (1)) is linear differential and, thus, it has several possible analytical solutions. One of them consists in applying Green’s Theorem to build a solution based on the distribution of sources and doublets on the problem boundaries (body and wake) yielding the following solution^{(Reference Katz and Plotkin29)}:

(2) $$\begin{equation} \Phi \left(x,y,z\right)=\frac{-1}{4\pi }\int _{S_B}\left[\delta \left(\frac{1}{r}\right)-\mu \; \textrm {n}\cdot \nabla \left(\frac{1}{r}\right)\right]dS + \Phi _{\infty }, \end{equation}$$

where δ and μ are, respectively, the source and doublet strength, r is the distance between the point (x, y, z) and the source/doublet location, $\textrm {n}$ is the normal vector to the surface, S_B the surface boundary and Φ_∞ the free-stream potential.

The numerical panel method based on the formulation described in Katz and Plotkin^{(Reference Katz and Plotkin29)} uses elements with constant intensity sources and doublets distributions. The application of these elements fulfills the far field boundary condition. However, to satisfy the impermeability condition, one has to impose one of three boundary conditions: (1) Neumann boundary condition; (2) Dirichlet boundary condition; or (3) mixed boundary condition. For the developed panel code, the Dirichlet boundary condition was used, which prescribes the total potential inside the surface as Φ = Φ_∞ and is the equivalent to state that the inner perturbation potential is zero. Equation (2) can then be expressed in the form:

(3) $$\begin{equation} \frac{1}{4\pi }\int _{body+wake}\mu \frac{\partial }{\partial \textrm {n}}\left(\frac{1}{r}\right)dS -\frac{1}{4\pi }\int _{body}\delta \left(\frac{1}{r}\right)dS = 0. \end{equation}$$

For this equation (Equation (3)) to have one particular solution, the source strength for solid bodies boundaries is prescribed to be:

(4) $$\begin{equation} \delta =-{\rm n}\cdot \nabla \Phi _{\infty }. \end{equation}$$

The wake is modelled as a constant doublet distribution, which strength is obtained by applying the Kutta condition:

(5) $$\begin{equation} \mu =\mu _u-\mu _l, \end{equation}$$

where μ _u and μ _l are, respectively, the doublet strengths at the upper and lower surfaces of trailing edge from which the wake element is shedding.

To solve this problem, one has to discretise the body surface and wake into several panels. The so-called aerodynamic influence coefficients were evaluated by computing the surface integrals in Equation (3) for each collocation point (placed at the centre of each panel). Since the source’s and doublets’ strengths are considered constant for each panel, the integrals depend only on the geometry. The doublet distribution (μ _k ) is differentiated using a second order central differentiation scheme to compute the tangential velocity (V) at the collocation points. The pressure (p) on each panel can be computed from Bernoulli’s equation for potential flow, assuming, as reference, the free-stream velocity (V _∞):

(6) $$\begin{equation} \frac{p_{\infty }-p}{\rho }=\frac{V^2}{2}-\frac{V_{\infty }^2}{2}+\frac{\partial \Phi }{\partial t}, \end{equation}$$

where p _∞ and ρ denote free-stream reference pressure and air density, respectively. The incompressible pressure coefficient (C _{p, 0}) can then be calculated for the steady case (∂Φ/∂t = 0) by,

(7) $$\begin{equation} C_{p,0} = \frac{p-p_{\infty }}{\frac{1}{2}\rho V_{\infty }^2} = 1 - \frac{V^2}{V_{\infty }^2}. \end{equation}$$

The pressure values for each panel are integrated over the surface to compute the total body pressure forces. Aiming at improving the accuracy of pressure estimation for higher subsonic flows, the Prandtl-Glauert rule^{(Reference Bertin11)} is applied, resulting in a corrected pressure coefficient,

(8) $$\begin{equation} C_{p} = \frac{C_{p,0}}{\sqrt{1-M_{\infty }^2}}, \end{equation}$$

where M _∞ denotes the free-stream Mach number. A friction coefficient (C_f ) and form factor (FF) for drag correction (C _Df ), considering the wing wetted and reference areas (S_wet and S_ref , respectively), was also implemented to account for viscous effects,

(9) $$\begin{equation} C_{D_f} = C_fFF\frac{S_{wet}}{S_{ref}}. \end{equation}$$

A detailed description of this correction (flat plate zero pressure gradient friction drag-based calculation) can be found in the literature^{(Reference Mason34)}.

2.1.3 Unsteady formulation

Despite the previous described methodology not accounting for time dependent terms, it is possible to perform small modifications to include them^{(Reference Katz and Plotkin29)}. Laplace’s equation (Equation (1)) is still valid for the unsteady case^{(Reference Katz and Plotkin29)}, although slight alterations are performed to its boundary conditions:

• • the impermeability condition is now called “zero normal flow on a solid surface” boundary condition and it depends on the local surface velocity;

• • the far field condition is maintained;

• • and, a new condition, known as Kelvin condition, depending on vortex strength, (which can be related to doublet strength) is required to determine the stream-wise variation with time of the strength of the vorticity shed into the wake.

For convenience, the unsteady motion is described in a body-fixed coordinate system on which “zero normal flow on a solid surface” boundary condition is applied. The body-fixed coordinate system is assumed to be known and is written in reference to the inertial frame. At t = 0, both coordinate systems coincide; while for t > 0, the relative motion of the body-fixed coordinate system origin is described by its location and orientation. The “zero normal flow on a solid surface” boundary condition can be rewritten in terms of the velocity (V ₀) and rate of rotation (Ω) of the body-fixed coordinate system:

(10) $$\begin{equation} \left(\nabla \Phi -V_0-v_{rel}-\Omega \times r\right)\cdot \textrm {n}=0, \end{equation}$$

with v_rel representing additional relative motion (e.g. in case of small oscillations or deformation speeds).

As for the steady case, the numerical panel methods use elements with constant intensity sources and doublets distributions, which fulfills the far field boundary condition. However, to fulfill the “zero normal flow on a solid surface,” the Dirichlet boundary condition can be applied as in the steady case. One can also rewrite the solution of Laplace’s equation similarly to Equation (3), but now the source strength for solid bodies becomes:

(11) $$\begin{equation} \delta =-{\rm n}\cdot \left(V_0+v_{rel}+\Omega \times r\right). \end{equation}$$

The wake model is not fixed as in the steady case; in fact, it increases in length (number of wake panels) with time and respects the following conditions:

• • no vorticity is added in the wake, Kelvin condition;

• • at each time step, the doublet strength is computed by applying Kutta condition (Equation (5)) at aerofoil trailing edges;

• • each wake panel length is computed by multiplying the time step (Δt) by the free-stream velocity (V _∞).

Considering now the time derivative in the Bernoulli’s equation (Equation (6)), one can calculate the incompressible pressure coefficient for the unsteady case,

(12) $$\begin{equation} C_{p,0} = 1 - \frac{V^2}{V_{\infty }^2}-\frac{2}{V_{\infty }^2}\frac{\partial \Phi }{\partial t}. \end{equation}$$

To solve this equation, a time discretisation method is required and since the objective is to compute the pressure coefficient at the time instance t + Δt, this method has to be implicit. One of the most commonly used, and the one implemented, is the Backward Euler Method^{(Reference Hirsch25)}, which is first order accurate and when applied to Equation (12) yields,

(13) $$\begin{equation} C_{p,0}^{t+\Delta t} = 1 - \frac{\left(V^{t+\Delta t}\right)^2}{V_{\infty }^2} - \frac{2}{V_{\infty }^2}\left(\frac{\Phi ^{t+\Delta t}-\Phi ^t}{\Delta t}\right). \end{equation}$$

The pressure values for each panel are integrated over the surface to compute the total body pressure forces. The compressibility and viscous drag corrections (Equation (8) and Equation (9), respectively) are also used in the unsteady formulation in order to improve the pressure and drag estimation accuracy.

2.1.4 Benchmarking

The implemented aerodynamic model (steady formulation) was compared in terms of chord-wise pressure coefficient, C_p , not only with experimental data^{(Reference Gregory and O’Reilly23)}, but also with numerical results, for a NACA 0012 aerofoil⁽³⁷⁾. Two turbulence models were used for this comparison: Spalart-Allmaras (SA)^{(Reference Spalart and Allmaras52)}; and standard Menter shear stress transport (SST)^{(Reference Menter35)}. For extracting the data in the developed framework, a wing with a high aspect-ratio (20 m of span and a chord of 1 m) was defined and the chord-wise C_P distribution was computed near the wing root for three different angles-of-attack: 0°, 10° and 15°. A Mach number of 0.2 (incompressible flow) and sea level altitude conditions were considered.

One can observe from Fig. 2 that the results obtained with the developed aerodynamic model are in good agreement with the experimental and numerical results for the chord-wise pressure coefficient distribution. A note should be made regarding the drag coefficient. As mentioned earlier, the Panel Method is inviscid, which leads to a lower drag coefficient prediction. It can be seen that for higher angles of attack (Fig. 2(c)) there is a slight deviation from the obtained chord-wise C_P distribution and the experimental values. This difference is expected to be more evident if C_D was plotted instead of C_P since the panel code does not properly take into account viscous effects, and at this angle-of-attack, there is flow separation.

Figure 2.

Comparison of the chord-wise pressure coefficient distribution with experimental^{(Reference Gregory and O’Reilly23)} and numerical results⁽³⁷⁾.

2.2.1 Problem formulation

The structural model is based on the finite element method (FEM) to describe the condensed 3D beam model. The implemented finite element is a 3D thin/slender beam element^{(Reference Soares50)} based on the Euler-Bernoulli beam theory^{(Reference Reddy45)} and composed of two nodes with six degrees-of-freedom (DOFs) each. This model consists of assembling four decoupled elements: one bar; two beams (one for each bending plane); and one torsion. The displacement field, {u}, which includes both displacements and rotations, can be defined for the assembled 3D beam element as a function of the shape function matrix, [N], and the nodal displacements, {d}, by the following expression:

(14) $$\begin{equation} \lbrace u\rbrace = [N] \lbrace d\rbrace . \end{equation}$$

The strain field, for this element, is composed of the axial strain, two bending curvatures (planes xy and xz) and twist angle; and can be defined in terms of strain-displacement matrix, [B], and nodal displacements as follows:

(15) $$\begin{equation} \lbrace \epsilon \rbrace = [B] \lbrace d\rbrace . \end{equation}$$

For this element, the stress field is also modified. It includes axial force and bending moments in the three planes. Furthermore, the material is assumed to be isotropic and elastic, which results in the following stress field:

(16) $$\begin{equation} \lbrace \sigma \rbrace = [D] \lbrace \epsilon \rbrace , \end{equation}$$

where [D] is the elasticity matrix. A brief description of the finite element method implemented in the structural model, whose formulation is based on the variational method, is presented next.

2.2.2 Static analysis

The potential energy, U, of the beam finite element is given by,

(17) $$\begin{equation} U = \frac{1}{2}\lbrace d\rbrace ^T[K]\lbrace d\rbrace - \lbrace d\rbrace ^T\lbrace F\rbrace , \end{equation}$$

where [K] and {F} are the stiffness matrix and total load vector of the element, respectively. The stiffness matrix and load vector are written in the element coordinate system. A transformation matrix is required to compute both this matrix and vector in the global coordinate system and allow assembling all the finite elements.

The static equilibrium is obtained by applying the principle of virtual work, yielding that the potential energy is stationary, with respect to the nodal displacements:

(18) $$\begin{equation} \frac{\partial U}{\partial \lbrace d\rbrace } = 0 \Rightarrow [K]\lbrace d\rbrace =\lbrace F\rbrace . \end{equation}$$

To solve this system of equations, an in-house linear solver that automatically selects the most suitable decomposition (LU, QR and Cholesky are possible) to invert the stiffness matrix is used.

2.2.3 Dynamic analysis

In a dynamic analysis, not only the elastic deformation is accounted for, but also the inertial properties are considered (the beam element model does not take into account damping). Therefore, it is necessary to compute the kinetic energy, T:

(19) $$\begin{equation} T = \frac{1}{2}{\lbrace \dot{d}\rbrace }^T[M]{\lbrace \dot{d}\rbrace }, \end{equation}$$

where [M] is the mass matrix of the element. The nodal velocity vector, ${\lbrace \dot{d}\rbrace }$ , is the time derivative of nodal displacement vector.

The potential energy is represented analogously to the one presented in the static analysis, although now the displacements and forces are in function of both time and space. For the dynamic equilibrium of the element in Lagrangian formulation, the equation is

(20) $$\begin{equation} \frac{d}{dt}\left(\frac{\partial T}{\partial {\lbrace \dot{d}\rbrace }}\right) + \frac{\partial U}{\partial \lbrace d\rbrace } = \left\lbrace 0\right\rbrace \Rightarrow [M]{\lbrace \ddot{d}\rbrace } + [K]\lbrace d\rbrace = \lbrace F\rbrace , \end{equation}$$

where the vector of accelerations, ${\lbrace \ddot{d}\rbrace }$ , is given by second order time derivative of the displacement vector. To solve this system, a time integration scheme is required and for the developed framework, both Newmark^{(Reference Bathe7)} and Bathe^{(Reference Bathe and Noh8)} methods were implemented. For the case of modal free-vibration analysis, there is no force applied; therefore, the equation of dynamic equilibrium is now an eigenvalue and eigenvector problem:

(21) $$\begin{equation} [K]\lbrace d\rbrace =\omega ^2[M]\lbrace d\rbrace . \end{equation}$$

An in-house eigensystem solver is used to compute the eigenvalues and eigenvectors of this problem.

2.2.4 Non-linear analyses

Aiming at enabling the model with the capability of incorporating geometric non-linearities, a few changes to both static and dynamic analyses were performed.

The stiffness matrix is now dependent of the deformed state; thus, in order to reach the non-linear static equilibrium at a given time step, t + Δt, an iterative method is required. The one implemented is the full Newton-Raphson method^{(Reference Bathe7)}, which aims at reaching the following equilibrium:

(22) $$\begin{equation} ^{t+\Delta t}\lbrace F\rbrace -^{t+\Delta t}\lbrace R\rbrace ^i=\lbrace 0\rbrace , \end{equation}$$

where the applied load, {F}, does not depend on the deformation, but the restoring force {R} does, with i being the current iteration. These restoring forces are the internal structural nodal forces, which are computed from the previous iteration stiffness matrix as follows:

(23) $$\begin{equation} ^{t+\Delta t}\lbrace R\rbrace ^i=^{t+\Delta t}[K]^{\left(i-1\right)}\lbrace d\rbrace ^i. \end{equation}$$

There are different approaches to treat large displacements and large rotations, although for the developed framework, NASTRAN™ methodology^{(Reference Lee30)} was followed. It consists in computing the net displacements and using the gimbal angle approach to treat, respectively, large displacements and large rotations. The implemented non-linear model does not account for large strains.

Since the main objective of this aeroelastic framework is to study high-aspect-ratio wings, where flexibility effects are of high importance, the element weight is also considered in the model. To do so, the gravity load, {F_g }, is computed by multiplying the mass matrix by the gravity acceleration vector, {g}:

(24) $$\begin{equation} \lbrace F_g\rbrace =[M]\lbrace g\rbrace . \end{equation}$$

Stress stiffening effects can also be modelled, if one considers a non-linear component on the strain-displacement matrix:

(25) $$\begin{equation} \lbrace \epsilon \rbrace = \left([B] + [B_{NL}]\right) \lbrace d\rbrace , \end{equation}$$

where the strain-displacement matrix has now a non-linear component, [B_NL ]. Besides the linear stiffness matrix, two more matrices can be computed: the large displacement stiffness matrix and the geometric stiffness matrix. The former is usually neglected^{(Reference Soares50)}, while the latter can be conveniently transformed as a function of stress state^{(Reference Zienkiewicz65)}.

2.2.5 Benchmarking

The structural model is a low-/medium-fidelity model consisting of an equivalent beam model, as already mentioned. For the purpose of validation, NASTRAN™ software was used. The problem devised for the benchmarking consists in a $10\,\text{m}$ aluminum (with a Young modulus of $70\,\text{GPa}$ , a density of $2,700\,\text{kg/m}^3$ and a shear modulus of $27\,\text{GPa}$ ) cantilever beam with a width and a height of $0.2\,\text{m}$ and $0.1\,\text{m}$ , respectively.

Different magnitudes of vertical loads were applied on the beam tip to evaluate the performance of the developed tool. The beam’s self-weight was accounted for in the calculations. Only constant forces were applied in the non-linear study. A uniform mesh of 500 co-linear CBEAM elements representing the cross-section was used for NASTRAN™ computations.

Both linear and non-linear static analyses presented very good agreement between the results from the MDO framework and NASTRAN™: a maximum error of 0.03% for the displacements and 0.26% for the maximum stress (for the load value) were reached for the non-linear solution. The main difference between NASTRAN™ and the implemented framework is in terms of non-linear solution convergence. While the developed tool updates the stiffness matrix in each iteration, using a full Newton-Raphson method^{(Reference Bathe7)}, NASTRAN™ does not perform a stiffness matrix update in every iteration, but employs a modified Newton-Raphson method instead. For this specific case, only one load step was required for both codes and both equilibrium and convergence were reached.

The same cantilever beam was also used in verifying the quality of the free-vibration analysis of the framework relative to the NASTRAN™ solution. Twelve mode shapes and the corresponding frequencies were extracted and compared: the relative errors for natural frequencies are very low (below 0.01%); and the mode shapes also exhibit very good agreement when one computes the coordinate modal assurance criterion (COMAC)^{(Reference Allemang3)} to compare both solutions (see Fig. 3).

Figure 3.

Comparison of the first five mode shapes using COMAC.

A structural dynamic test case was also devised to compare the performance of linear and non-linear dynamic solvers. This test case consists of a cantilever beam (with a length of 0.2540 m and a square cross-section with $0.0254\,\text{m}$ sides) with a uniform load (of 19, 700 Pa) applied along its length and constant in time. An isotropic material with a Young's modulus of 82.737 MPa, a density of $10.943\;\text{kg/m}^3$ and a Poison ratio of ν = 0.2 was considered. Both Newmark (trapezoidal rule) and Bathe methods were compared with results found in the literature using the same time step ( $\Delta t = 0.00015\,\text{s}$ ) and mesh (20 finite elements). Despite the beam shape, literature models are based on plates. So, in order to expand the comparison to other elements, linear and non-linear transient analyses were carried out in ANSYS™ using the same time step and a mesh of solid elements (SOLID187) with a 80 × 8 × 8 (length, width and height) discretisation. In Fig. 4, the different linear and non-linear transient responses measured at the beam tip are depicted for the first $0.006\,\text{s}$ .

Figure 4.

Structural dynamic response, literature results are marked as Relvas^{(Reference Relvas46)} and Belytschko et al^{(Reference Belytschko, Lin and Chen-Shyh10)} in the legends.

One can observe from Fig. 4 that both linear and non-linear framework solutions are in very good agreement with those found in the literature. Maximum errors under 2% and 3% were verified for linear and non-linear solutions, respectively. Despite this good agreement with the literature results, the non-linear dynamic solutions obtained using the developed framework are different from those of the solid element, mainly in terms of amplitude, which is lower (maximum errors of 8.07% and 8.26% were observed when comparing the framework results obtained using the Newmark and Bathe methods, respectively). It is worth noting that the implemented beam finite element has a different formulation (based on Euler-Bernoulli beam theory) to those of the solid element, which results in a more accentuated difference in the non-linear solution when very large deformations occur (the vertical tip displacement corresponds to more than 70% of the total beam length for the linear solution). The linear and non-linear vertical tip deflection values are summarised in Table 1.

Table 1

Vertical tip deflection

2.3.1 FSI

The FSI methodology implemented can be classified as: (1) monolithic^{(Reference Hulshoff27)}, since the solution advances with the convergence between structural and aerodynamic solutions at a time instant in the unsteady case (or at a given load step for the steady case); and (2) loosely coupled^{(Reference Bazilevs, Takizawa and Tezduyar9)} since the fluid and structure equations are solved separately and coupled in an interface model until convergence is reached.

This algorithm is responsible for transferring the loads from the aerodynamic panels to the structure and transferring the structural displacements to the aerodynamic panels. Usually, accurate modelling of the fluid dynamics requires a detailed mesh, while for modelling the structural dynamics, a relatively coarse mesh suffices^{(Reference Ballmann5)}. As a consequence, the two meshes are not coincidental and, therefore, a proper load and displacement transfer method is necessary.

The aerodynamic loads computed at each aerodynamic node are transferred to the structural mesh by means of generating an equivalent load system, which is applied to the closest beam finite element. Despite ensuring consistency, this approach does not guarantee energy conservation. Nevertheless, for flexible structures, such as high-aspect-ratio wings, if one ensures energy conservation, one might obtain large non-physical oscillations in the pressure force resulting in a significant negative influence on results accuracy^{(Reference Rixen47)}.

No aerodynamic mesh deformation is performed; instead, a new mesh of the deformed shape is generated for each aerodynamic calculation. This mesh is built based on the structural displacements, which are transferred to aerodynamic nodes by assuming the aerofoil cross-section as rigid, i.e. the structural nodal displacement is transmitted to the aerodynamic node by a rigid-link.

The FSI routine is an iterative process which converges when equilibrium is reached, i.e. when both displacements and loads have converged. In the steady case, the aerodynamic model is re-meshed at each equilibrium iteration. In order to facilitate convergence, a few features were included in the developed framework: firstly, in the case where convergence has not been reached, the FSI module performs a bisection of the load and displacement and starts at the previous equilibrium state (which, in the case of the bisection taking place in the first iteration, is the undeformed state); secondly, it is possible to begin with just a fraction of the total load meant to be transferred to the structure. The convergence criteria is given by,

(26) $$\begin{equation} \frac{1}{NN}\sum _{i=1}^{NN}\left|\frac{d_i^n-d_i^{n-1}}{d_i^n}\right|\le 1\times 10^{-7}, \end{equation}$$

where NN and d_i are, respectively, the number of structural nodes and nodal displacement (displacements and rotations) i, n is the current load step and n − 1 the previous load step. Figure 5 shows the implemented FSI flowchart for static aeroelastic calculations.

Figure 5.

Steady FSI algorithm flowchart.

All structural time marching schemes assume the loading as a function of time is known. For an unsteady FSI analysis, this is not the case, since the aerodynamic loading at a particular time is determined by the deformed position of the wing/aircraft and by the local deformation speed influencing, for instance, the local angle-of-attack. In particular, when using an unsteady panel method, such as the one implemented, the aerodynamic loading calculation requires a re-meshing or mesh deformation procedure according to the aircraft deformation and that the local velocity at the collocation points is known. Thus, an iterative process is done in order to match the load applied to the structure with the aerodynamic load computed at the end of the load step. Figure 6 shows the implemented time marching procedure; the iterative loop box is depicted in Fig. 7.

Figure 6.

Unsteady FSI scheme.

Figure 7.

Unsteady FSI iterative loop.

For an initial load guess (usually equal to the load at the end of the previous time step), the procedure uses a Newmark or Bathe time integration method in the structural solver for the calculation of the nodes displacements, velocities and accelerations at the end of the time step.

With the nodal displacements and displacement velocities available, one of two possible work flows follows, depending on the model being constrained or not. If the model is free, an intermediate step is performed in order to obtain information about the Centre of Gravity (CG) position, rigid body orientation and linear and rotation speeds by decomposing the total displacement and displacement velocities into rigid+deformation components, as depicted in Fig. 8. This decomposition allows the use of the framework for rigid body flight simulation by neglecting structural deformation. With the resultant loading from aerodynamic and other sources (e.g. propulsion and gravity), a 4th-order Runge-Kutta integration scheme with an embedded predictor-corrector scheme for aero-structural convergence is used to determine the rigid body CG position and orientation, and CG linear and rotation speeds at the end of the time step, followed by a predictor-corrector scheme used for aero-structural convergence and to determine the linear and rotational accelerations, as well as in the end of the time step. If the model is constrained (e.g. half wing clamped at the root, as is the case studied in this work) the orientation of the model is assumed fixed and the rigid body components are assumed null. In this case, the predictor-corrector scheme is used for aero-structural convergence before advancing the solution in time.

Figure 8.

Example of mesh deformation calculation for a CG position displacement from time t to time t = t + Δt.

2.3.2 Re-meshing procedure

The re-meshing procedure is based on rigid links between structural and aerodynamic nodes, as already mentioned. The aerodynamic mesh is constructed by defining a set of points and associated aerofoil (or fuselage) shapes and shape-related parameters (e.g. chord, twist and wing-box chord-wise limits), in a piece-wise linear or other shape function approximation. These points are the ones to which the structural displacement is passed on to by a rigid link connecting them to the associated structural nodes. In the present case, the aerofoil shape is assumed to be rigid and, therefore, it is translated and rotated according to the structural deformation calculated. In the unsteady case, not only the mesh is updated according to the structural deformation but also the deformation speed is calculated for each aerodynamic panel control point, by linear interpolation of the deformation speeds calculated for the construction points in which the aerodynamic panel is between. This relative velocity is then added to zero normal flow boundary condition for the aerodynamic panel.

2.3.3 Static aeroelastic analysis

For the static aeroelastic calculations, the steady FSI algorithm, shown in Fig. 5, is employed, which calculates the linear or non-linear static equilibrium position.

2.3.4 Dynamic aeroelastic analysis

Two dynamic aeroelastic methodologies were implemented: (1) frequency domain and (2) time domain. The former consists in calculating the flutter boundary at the linear or non-linear equilibrium state (computed with the steady FSI algorithm) using the p-k method^{(Reference Hodges and Pierce26)}. In the latter, the unsteady FSI (time-marching scheme, shown in Fig. 6) is used to compute the non-linear aeroelastic deformed state in time.

NASTRAN™ was used for flutter speed calculation in the frequency domain. In this methodology, only geometric non-linearities are transmitted to the NASTRAN™ model through the generation of a deformed model and the flutter speed is computed for the static aeroelastic non-linear equilibrium state.

2.3.5 Benchmarking

An experimental aeroelastic model of a high-aspect-ratio wing with a slender body, presented by Tang and Dowell^{(Reference Tang and Dowell57)}, has been used for benchmarking non-linear aeroelastic numerical codes and was, likewise, applied here to benchmark the implemented FSI algorithm: a rectangular spar with the same inertia properties was generated with the developed tool to represent the internal structure of the wing and the slender body was modelled as a mass point with equivalent mass and inertias. The generated aerodynamic model does not take into account the aerodynamics of the slender body.

Prior to performing aeroelastic analysis, modal identification took place and the results were compared with the literature and experimental results^{(Reference Arena, Lacarbonara and Marzocca4)}. Good agreement with the numerical data of Arena et al^{(Reference Arena, Lacarbonara and Marzocca4)} was verified, as well as with the experimental data, as can be observed in Table 2. That applies to the model with and without the slender body. It is clear that, although the addition of the slender body at the wing tip reduces all natural frequencies, the first torsion and the first chord bending modes are the ones most affected.

Table 2

First 5 natural frequencies for the Framework (Frame.), Experimental (Exp.) and Arena et al^{(Reference Arena, Lacarbonara and Marzocca4)} numerical results, with and without slender body. The Δ values represent the relative errors

In this test case, non-linear aeroelastic solutions were computed for ten flow speeds (ranging from 5 m/s to 35 m/s), at sea level altitude and for an angle-of-attack of 1°. The results are shown in Fig. 9.

Figure 9.

Comparison of vertical tip displacement and twist angle as a function of the flow speed with experimental and numerical results of Tang and Dowell^{(Reference Tang and Dowell57)}.

One can observe from Fig. 9(a) that the vertical tip displacement is close to the numerical and experimental data for flow speeds of up to 30 m/s. For higher speeds, the tip displacement is small, and the FSI module tends to converge faster, which together with the fact that the slender body was not included in the aerodynamic model leads to higher errors in the vertical tip displacement, when compared with the experimental data. Good agreement was found, despite the differences in the aerodynamic models (2D ONERA aerodynamic model^{(Reference Petot43)}). It is noticeable the increase of the difference between experimental and numerical tip twist angle, i.e. when the aerodynamic of the slender body becomes more important which causes an increase of the effective angle-of-attack and, thus, an increase in vertical tip displacement. The absence of the slender body in the aerodynamic model explains the discrepancies between the developed framework results and the experimental data. In fact, the structural model closely represents the structure tested in the wind tunnel, although the absence of the slender body in the aerodynamic prevent a direct comparison between the results.