Aeroelastic simulations in OpenVOGEL

By coupling the unsteady vortex lattice method to a simple finite elements-based elastic model, OpenVOGEL tries to provide a basic approach to the solution of several complex problems associated to elastic lifting surfaces, namely:

change on aerodynamic properties due to wing deformation (change on lift slope and moment slope)
divergence
aileron reversal
response to gusts
flutter (resonance of aeroelastic nature)

These problems are important to be understood since they can bring the aircraft to unsafe situations with even catastrophic consequences. The first three problems belong to the field of "static aeroelasticity", while the last two belong to the "dynamic aeroelasticity". Note that these last two problems are not yet fine tuned in our calculation model.

Representation of a typical aero-servo-elastic problem: the reduction of the lift coefficient with airspeed due to twisting in a flexible clamped-free wing with a flap down. The velocity at which the lift is no longer increased when deflecting the flap is called *aileron reversal airspeed*.

Fluttering wing
Wing in divergent flutter.

What is aeroelasticity, and how the simulation works

Aeroelasticity is a very wide topic, which basically treats the interaction between a fluid (air) and an elastic boundary, that is to say, a boundary that is able to deform when subjected to loads. As the fluid flows it loads the boundary, as consequence the boundary deforms, and in return the flow is modified. The whole process is in permanent feedback, and therefore we talk about an interaction. VOGEL is only oriented to the static and dynamic response of aircraft wings due to air loads. If this topic is new for you, I would recommend you to take a look at the following book before even starting: Introduction to Structural Dynamics and Aeroelasticity, Hodges & Pierce, Cambridge. As you will note there, solving aeroelastic problems on wings is not always an easy task. First of all, one need to count with an algorithm based on a mathematical model capable of providing the unsteady air loads (steady aerodynamics is most of the times not suitable to solve this kind of problems). Secondly, the aircraft structure has to be modeled somehow, and finally a link between the motion of the structure and the airloads has to be established conforming to the laws of motion. By using the UVLM as the unsteady-aerodynamics solver platform, the only "extras" in VOGEL are the solution to the structural dynamic problem and the introduction of an effective link coupling both problems.

Structural model and aeroelastic link

There are many options to model the structure of the wing. The most practical way to do it is by finite element methods. There are several options:

working with 3D solid elements.
working with 3D shell elements (DKT elements, Reissner-Midlin or any other plate theory).
working with simple 3D beam elements.
working with a mix of elements.

Since the aerodynamic problem already takes many degrees of freedom, the third option can be very efficient, and therefore it is the one that has been selected. VOGEL uses a mesh of simple 2-noded 3D beam elements to model the wing structure. By doing this a reduced number of degrees of freedom is obtained, and a simple but very effective link can be implemented. The elastic and inertial properties of small portions of the wing are thus concentrated on associated beam elements distributed span-wise over the wing.

Besides the element type, another very important thing is the method to be used to solve the dynamic problem (the equations of motion). There are a couple of options to do that too:

Direct time integration method
Modal decomposition method.

The second method is the one implemented in VOGEL. Although it requires some initial effort to calculate the dynamic modes, once that problem is solved, it is much easier to handle and it requires less memory usage. The main disadvantage is that It is only valid under linear conditions and low deformations. This means that if we are for instance planning to go behind the skin buckling deformation, it will not provide accurate results anymore. However, for most of the problems it does provide great results.

The aeroelastic link

An "aeroelastic link" is a sort of system coupling on one side the aerodynamic loads on the lattice to the structure and, on the other side, the structure response to the lattice motion.

The aeroelastic linking consists of three basic steps:

The air-loads from certain rings of the vortex lattice are sent to certain nodes on the structure.
As result of these actions (a force and a moment on each structural node), and according to Newton's laws of motion, the structure undergoes motion.
The motions of the structural degrees of freedom modifies the geometry of the vortex lattice, and the loop is repeated.

It is natural to think that this can best be done by splitting the problem in two parts by means of the following object types:

Mechanic links: each holding a group of chord-wise ring-stripes from the vortex lattice associated to a beam element (refer to the figure below to get a better idea).
Kinematic links: each holding one structural node linked to a group of chord-wise nodal points on the vortex lattice. We will assume that these kinematic links act in a rigid way, that is to say, that every displacement and rotation on the structural node is sent to the lattice linked nodes as if they were rigidly connected.

All the required information can be held in sets of the above object types. In order to bring this idea into practice, VOGEL handles them in two apart .NET classes holding reference to a bunch of data associated to the lattices and the structure. Both classes are loaded prior calculation with the required elements, and later on used to update the forces and the motion. This is explained next.

Mechanic links

The MechaLink class holds one structural beam element and a list of vortex rings. It contains a very important method: "TransferLoads()". This method uses the location of each control point and the aerodynamic load on each linked vortex ring to build a total load vector which is assumed to be equally distributed between both extremes of the beam element.

Kinematic links

The KinematicLink class holds a list of chord wise lattice nodal points and only one structural nodal point. It contains the method "TransferMotion()" which, based on the structural displacement and velocity of the structural node, updates the position of the lattice points. To do this, this class holds the initial position of the lattice points in an apart list, so that displacements can be added to the original position every time an update is required.

As I have said before, lattice nodes will follow the displacement of their associated structural node according to a simple kinematic law: a "rigid link". A different law could be considered, although it does not make any sense at this point to make the problem even more difficult without measurements showing how it better should be.

Time integration scheme

By introducing the required MechaLinks and KinematicLinks, a great deal of the problem has been properly solved. These two classes will provide all required information for the aeroelastic coupling in a very well organized manner. The most important part of the problem is however still pending, that is the algorithm through which the equilibrium states will be found. There are here several ways to proceed, and choosing one or another will depend on what kind of analysis we are interested in.

If we are are only interested in the static steady state, then we can choose a simple explicit time integration. Algorithms belonging to this kind never use the new predicted state to feed themselves back in other to converge to a new equilibrium state. Algorithm 1 shown below belongs to this kind. If the aeroelastic transit tends to a static steady state, this algorithm will guarantee equilibrium there and will be very effective.

If we are concerned about the transit states, then we need a more sophisticated algorithm, because we need to seek for equilibrium at each time step. In this case we need an implicit algorithm that recalculates the air-loads and displacements several times for each time step until the equations of motion agree, that is to say, until the loads in the predicted position are in dynamic equilibrium with the predicted motion after the considered time step. One way to do that is by following a Gauss-Seidel scheme. Algorithms 2 and 3 are of this kind.

Algorithm 1: Explicit time integration

Calculate RHS
Build matrix
Calculate circulation
Calculate airloads
Transfer loads to structure
Predict modal motion with new load
Transfer motion to lattice
Calculate velocity on wake nodal points
Update wakes

Algorithm 2: Gauss-Seidel scheme with implicit time integration and preconditioned wakes

If t > 1
1. Calculate velocity on wake nodal points
2. Update wakes
Find equilibrium state at t+1:
1. Calculate RHS
2. Build matrix
3. Calculate circulation
4. Calculate air-loads
5. Transfer loads to structure
6. Predict modal motion with new load
7. Transfer motion to lattice
8. If not converged go back to a.

Algorithm 3: Gauss-Seidel scheme with implicit time integration and updated wakes (current method)

Find equilibrium state at t+1
1. If it is not the first implicit step, reestablish wake to previous time step.
2. Calculate velocity on wake nodal points
3. Update wakes
4. Calculate RHS
5. Build matrix
6. Calculate circulation
7. Calculate airloads
8. Transfer loads to structure
9. Predict modal motion with new load
10. Transfer motion to lattice
11. If not converged go back to a.

Method for the integration of the equations of motion

The integration of the uncoupled equation of motion requires of a proper algorithm. There are many algorithms that can be useful here, as described in "Finite Element Procedures" (Bathe, 2006). Due to its simplicity, I started working with the "central difference method". However, that and others methods are not always stable, and therefore, convergence problems might be found for certain dynamic modes under certain time-steps. Newmark method is the most suitable for this problem because it proves to be unconditionally stable for any time-step, and because it doesn't introduce numerical damping (although it does shrink or expand the frequency for too-large time steps).

The problem of matching time steps

A very important fact on the integration of the equations of motion is that they are solved simultaneously with the UVLM, and therefore time steps of both methods should somehow match. In general, the integration of the structural equations will require a different time step from that one on the UVLM, and therefore a partition of one of them will be required. The longest required time step can therefore be partitioned into equally spaced steps, and this refinement can be used to integrate the other set of equations. Say for instance that the UVLM requires of a 0.01s time step in order to get a regular mesh on the wakes (this condition being dependent on the length of the chord-wise vortices and the airspeed) , and that the equation associated to the highest dynamic mode requires of a time step of 0.001s in order to get an accurate response. In that case, we could run exactly 10 structural time steps per aerodynamic step. For most accuracy, the aerodynamic load should be updated on each structural time-step, however this requires an huge calculation effort. Moreover, it would not guarantee more accuracy if the circulation on the wakes is not properly updated. Our way to solve the problem is thus suitable for when there is a soft aeroelastic coupling. There are other advanced techniques to achieve a good aeroelastic link, however, they all require iteratively updating the aerodynamic loads, something very costly from the computational point of view.