OpenVOGEL/Printable version

OpenVOGEL

The current, editable version of this book is available in Wikibooks, the open-content textbooks collection, at
https://en.wikibooks.org/wiki/OpenVOGEL

Permission is granted to copy, distribute, and/or modify this document under the terms of the Creative Commons Attribution-ShareAlike 3.0 License.

Introduction

Foreword

OpenVOGEL is an open source project founded with as goal to provide free access to a computer program that would allow the numerical study of aeromechanic problems (aerodynamics + elasticity + dynamics). OpenVOGEL can be used to create from scratch, calculate and analyse several aspects of an aircraft model. The software integrates grid generators, unsteady flow theory based in first order panels, structural dynamics by finite elements (modal decomposition) and a graphical user interface.

OpenVOGEL relies in a series of common software packages that are implemented in two separate user applications: Tucan (a user friendly GUI) and the Console (a command line tool). Throughout this Wikibook you will find information about what these two programs are capable of, how they works, and how they can be used in aircraft design.

Aircraft model built and simulated with OpenVOGEL Tucan

OpenVOGEL is open source, which means that everyone has access to the source code. The source code is published under the General Public License (GPLv3). To help with the development of the code, please visit repository.

Note to newcomers

If you love airplanes but you are new to the world of aeromechanics, take into account that the topic is difficult in many ways. If you don't have a solid background in basic math (algebra and analysis), you might have a very difficult time trying to figure out what this is all about. I do not consider myself an expert in mathematics so you don't have to be one either, but the basics like working with Euclidean vectors you cannot miss. To understand how an aircraft flies you also need a basic knowledge of fluid dynamics and, in general, classic (Newtonian) mechanics.

This project targets mainly aerospace engineers and undergraduate students actively working on aircraft design. However, this does not prevent amateurs or less qualified people from using the program as guidance to understand how to fly or design an RC model or drone. I have tried to make the program accessible to a broad public.

Over this wikibook

The intention of this wikibook is to make the documentation of the project as rich as possible through the contribution of users. This wikibook aims to cover theoretic notions behind the software and information about how to use it, how NOT to use it, and how to develop it. Open source code is great, but to get the most of it the code must be clearly specified and documented. Many open source projects have become incredibly popular, but their lack of documentation causes a lot of confusion and wrong implementations. To avoid this, we try to provide here as much information as possible so that users can have a good idea of both the bigger picture and the low level math.

This book is currently divided in the following chapters:

• User's guide
• Source code
• Aerodynamics
• Aeroelasticity
• Free flight
• Validation

The evolution of the program

Sometimes it is useful to know a little bit of history to understand why something is the way it is, so I would like to write here a brief introduction about the evolution of OpenVOGEL.

OpenVOGEL was born at the end of the year 2009 as a small research project during my last years at university. It was first developed in FORTRAN, and it was actually only intended to study the behavior of droplets being sprayed by fumigation airplanes flying at very low altitude, a curiosity that can only emerge in a country where corn and soja fields cover an extent of gigantic proportions.

As the program matured, I started to develop a Human-Machine Interface (HMI) for it in Visual Basic 6; unsatisfied with the result and having learned VB.NET at work, I decided to rewrite the whole calculation core in a full object-oriented fashion using Visual Studio and .NET. This compiler (actually virtual machine) was chosen because it provides access to innumerable useful standard libraries and because it is stable and mature enough for doing scientific calculations. One of the advantages is that it saves lot of development time, mainly because memory is managed by something called a garbage collector, which is not present in languages like C, Ada or Fortran. This makes life easier for the developer by letting them focus on the functionality of the software instead of in how the computer should do the job. This is not always necessary to have and for some applications such as real-time applications, it might even be undesired. When it comes about prototyping science projects, however, it clearly is an advantage. Of course this is paid with some loss of performance but for the purpose of this project, the balance was clearly in favor of easier coding. Within .NET there are many compliant languages, Visual Basic (VB) and C# being the most common two. I chose VB instead if C# simply because it is more natural to read and easier to learn. However, many libraries are found in C# and have been linked to the project without problems.

When it comes to coding, the evolution of OpenVOGEL in .NET is naturally linked to the evolution of Visual Studio Community Edition. Since 2010, Microsoft has been publishing a free version of Visual Studio which has evolved rapidly and in favor of open source communities. It started with the Express editions and evolved towards a series of Community editions, which are a lot more feature-rich than the former but with extra license limitations. In 2014 Microsoft pushed .NET as a standard specification, and since then it has been implemented in a parallel open source project called Mono. This project also includes an open source editor called Mono Develop which would eventually be used to develop OpenVOGEL calculation core and make it cross platform.

Let's now talk a bit more about the .NET versions of the program. A program is not made from one day to another. Normally, the life cycle is as follows. The program starts with something simple, evolves to something complex, and then reaches a point where it becomes so complex that it has to be restructured from scratch again. Only when you know exactly from the beginning what the program should be capable of are you in a position to make a strong and flexible kernel. With OpenVOGEL it happened just like that. The current version is the result of rewriting parts to make them more general. The first version of the program only included slender lifting surfaces, and it was only able of handling Neumann boundary conditions (by explicitly fixing zero normal velocity at the boundary). After validation of the results, the aeroelastic algorithms followed. The development of Dirichlet boundary conditions (zero internal velocity potential) took more time to develop, and it is one of the latest implemented features. At least the first four years of development were done in private because the program was not yet mature enough to meet the world. The program first saw light when it became open source in the year 2015 with the idea of making science instead of marketing, and because I sensed the lack of a free friendly program in this area. While the market of proprietary software is very rich, there is very little up-to-date open source software based in panel methods. In fact, most books still present popular FORTRAN routines of the 70's, when object-oriented programming was still in the early development phase. Based on today's technology, OpenVOGEL offers a new perspective on aeromechanical modeling.

Limitations and capabilities

As we will see in the chapters of this book, the aerodynamic method that is implemented in the software can be very powerful for solving some problems and insufficient for others. The restrictions are of differing natures: due to complex geometry, due to the simplified representation of true physics, or even due to limited calculation power. The restrictions of the software will be pointed out as they become evident while describing the different modules of the software. However, I would like to make a small resume of the most important facts. Let me start with a brief description of what simulation means.

Simulating a system is trying to predict its behavior by watching the behavior of a model, which can be a material object or a collection of linked algorithms. Simulating a system with computer algorithms is clearly not like a real test, where real fluids are used and the real equations are naturally integrated. For the last three decades or more, computer simulations have improved dramatically but still face limitations. The problem of numerical simulations is treating the complex math that represents the system behavior with a limited amount of computation resources. We cannot keep track of all molecules in a fluid because there is no sufficient storage in a computer to do so, so the problem needs to be simplified and condensed in a process called discretization.

The problem of real simulations, on the other hand, is that the environment has to be carefully adjusted and controlled to match the target situation. Even if using the same fluid, a real model can deviate from the target situation. For example, consider testing an airplane using a wind tunnel. First you need a model of the plane matching the shape and roughness of the real model. Then you need the wind tunnel, which due to the limited volume does not behave the same as an open fluid domain. Boundary layers develop in the walls, and they introduce effects that are not encountered in free flight. Moreover, the airspeed and the static pressure need to be adjusted to match the Reynolds number, and this could be a challenge. Also, if the flow is not well controlled downstream, at the entrance of the tunnel, turbulence could be introduced that is unreal or out of proportion. Finally, you will need an accurate instrument to measure the air loads. And this is only to mention a few of the many possible sources of deviation.

So we could conclude that simulating a system is never exactly like reality, and no matter the method that is used it is always a challenge. It is an approach that could accurately represent some aspects of the target situation, while at the same time failing to predict others.

Probably the most important thing to make clear in our case is that OpenVOGEL is totally based on potential flow theory. This is an idealization made in fluid dynamics. It does not match reality, but it simplifies the mathematical description of the problem. For many reasons, this simplification must be done.

The potential flow theory states that when viscosity is zero (or, in practice, when the Reynolds number tends to infinity), the flow becomes irrotational and subsequently all we need to describe the motion is a potential function in the form of an unsteady scalar field. The velocity becomes the gradient of that field, and the pressure is linked to the square of the velocity by the very well known Bernoulli's equation. It is a beautiful and understandable theory for sure. However, there are no skin friction forces in potential flow theory, not to mention micro turbulence. So you will never get an accurate description of how boundary layers behave. You will never even get to see one. I know this sounds very disappointing, but don't quit reading yet because in practice it manages to give very good predictions of the main aerodynamic forces and in just a couple of minutes.

Other CFD programs are able to solve Navier Stokes equations, and this brings some advantages. For instance, these programs usually include friction forces and turbulence models. However, the algorithms behind these programs are far more complex since they are typically based in finite elements or finite difference methods. They require many more degrees of freedom and a more complex geometrical description, and therefore they demand more memory and CPU capacity. They also need to be fine tuned to work, usually by modifying control parameters that are adjusted to match real simulations. Because of this, most designers will not use them in the early design of an aircraft where different geometry configurations need to be quickly compared. Complex CFD programs are more natural to be used in research of cutting edge technologies, where accuracy plays an important role. So you are probably beginning to realize from these thoughts the importance of using simple methods, such as the one we will study here, for common engineering applications.

User's guide

The following tutorial pages are intended for end users. Hence, here you will not find information about how the program has been built, but about how it works and how to use it in practical aircraft design.

Graphic interface

Part 1: simulating a simple wing model

In this guide we show you how to:

• Create and edit a lifting surface.
• Load polar curves to simulate skin drag using XFoil or other data sources.
• Start a steady simulation.

Part 2: creating a fuselage

In this guide we show you how to:

• Create and edit a native OpenVOGEL fuselage.
• Import fuselages from OpenVSP or other sources.

Part 3: propellers

In this guide we show you how to:

• Create and edit a propeller.

Calculation console

Part 1: introduction to the Console

In this guide we show you how to:

• Introduce commands and work with the cosole.
• Run the simulation remotely from Tucan.
• Setup Intel MKL for improved calculation performance.

Part 2: running the simulation from the console

In this guide we show you how to:

• Run a steady analysis in the Console.
• Run different standard batch simulations and interpret the results.

SourceCode

Source code

Introduction

As it was said in the introduction, OpenVOGEL has been programmed in the .NET frameworks, mostly using Visual Basics. External libraries written in C# have been linked, but their development is out of the scope of this project, since they have remained almost untouched. I am conscious that it is not always easy to understand the source code written by someone else, so in this chapter I will try to make an effort to explain in words and as clear as possible how it all has been organized. If you had experience in FORTRAN code (such as PanAir) and you want to learn OpenVOGEL, take into account that object oriented programming consists in a very different approach. Hopefully, after reading this you will be able to find the way through the code and make your own adaptations the way you like.

Libraries

In OpenVOGEL the code is distributed in different libraries, and there is a logic for this. To understand the structure of libraries and their relationship, follow the diagram here next.

Architecture at large

All low level mathematical procedures are located in OpenVOGEL.DotNumerics. This is a fork of the DotNumerics project created by Jose Antonio De Santiago Castillo, which is basically an automatic translation of LAPACK and BLAS from FORTRAN to C#. In this project I have added the Subspace Iteration method (Bathe), which is a very specific algorithm suitable for finding the lowest eigen-values and vectors in systems of the type Mv=aKv with a large number of degrees of freedom. Additionally, I have added to the library bindings to the Intel MKL that can be used optionally for improved calculation performance instead of the native .NET routines. These bindings retain however the DotNumerics high level API. DotNumerics is thus basically used for solving the system of linear equations by LU decomposition, and to find the vibration modes of the wings.

OpenVOGEL has also its own linear algebra package, which is stored in the library called OpenVOGEL.Math. Here I have introduced vectors, numeric integration and other useful objects that are essential for the aerodynamic, dynamic and structural algorithms. To have an idea of how important this is, think that all vectors in the project are either Vector2 or Vector3 classes residing in OpenVOGEL.Math.EuclideanSpace.

The actual potential flow and aeroelastic solver is contained in the OpenVOGEL.AeroTools library. The calculation core only uses the previously mentioned libraries, and contains general definitions that can be embedded in any external project without the necessity of the Tucan or the Console.

In addition to calculation libraries the project also offers modeling tools that can generate the mesh of airplane models through a parametric description of the geometry. This is contained in the OpenVOGEL.DesignTools library. This library also offers the conversion procedures that are necessary to go from the design model to the calculation model and vice versa. Again, this library can be embedded in any external project without requiring Tucan or the Console.

The OpenVOGEL.Console project is a console application that can easily run in Windows or Linux (using Mono). This console is able to read, calculate and write any kind of OpenVOGEL project, just as Tucan does it (i.e. calling exactly the same procedures), but without graphical interface. This application is offered for advanced users that are able to do the pre- and post-processing by themselves. You can use it, for instance, to run a batch analysis and to perform your own customized analysis.

Finally, OpenVOGEL.Tucan is our main solution targeted for the general public. It relies on System.Drawing, Winforms and OpenGL to generate a graphical representation of the data. This solution is currently only available for Windows (7 and 10).

Classes

Now you understand the architecture of the project, it is time to take a look at the source code. To understand the source code, you first need to know about object orientation and how .NET implements it, since most of the code has been written that way. If you are new to this, then I would recommend you start by reading some dedicated book (there are many of them), along with the documentation of .NET provided by Microsoft, which is very rich.

In short, the difference between object orientation and procedural code is that in the former style we focus more our attention in how the system is divided in functional components, and how these components work together as an assembly. In procedural code, on the other hand, you would focus more on the different actions that a system performs to reach a goal. The data is then passed from one routine to another in place of residing inside an object.

In object oriented programming you would dismantle your system into several types of small components called classes, and then instantiate these classes in as many components as you need to build your assembly. Each class encapsulates data in the form of fields and properties, and is able to perform actions through functions and procedures, also called methods (try to remember all these concepts, since they constitute the soul of object orientation). Moreover, components that share some properties or procedures can be built to inherit a common ancestor, something that is called inheritance. For more clarity, let me put this idea in a down-to-earth example. In OpenVOGEL we need to model airplanes, which are basically made of several wings, fuselajes and nacelles. However, it is clear that all these three components are in fact some kind of surface, so they share in fact a set of common properties. They all have a mesh, and they all can be selected, moved, rotated and written into a file. So independently on what kind of surface they are, they all inherit a primitive class called Surface.

Public MustInherit Class Surface 
   Implements IOperational 
   Implements ISelectable
   [...]
End Class

As you can see, Surface is forced to implement the interfaces IOperational, and ISelectable, which will guarantee that all surfaces will expose a set of basic methods for handling them. The interface IOperational will force the overriding class to implement rotation and translation routines, and the ISelectable interface will force the overriding class to implement 3D selection routines. So you see that by implementing these interfaces we are declaring a consistent way of working.

Among the advantages of this way of thinking is clarity in code reusability. For example, lets imagine that we need a new kind of lifting surface, lets say, a VerticalEmpennage. If we would not have classified all our surfaces under a common ancestor, then Vertical empennage should be generated from scratch, without profiting from the general procedures of a normal LiftingSurface, or a Surface. But by letting all surfaces be a kind of Surface, then we can threat it exactly the same way as all others. VerticalEmpennage could then be derived from LiftingSurface, implement all its properties and methods, but then provide an interface more oriented towards the properties of an empennage, like for instance the configuration of the rudder.

Basic coding principles

All these tools provided by the object oriented technique are in fact meant to make your code work as a machine. If your interest in life diverges to mechanics, then maybe programming is not your forte. So to get the things right from the beginning, there is a fundamental principle of coding that you should keep in mind. Programming is not very different from building a hardware machine. A successful machine you build by assembling units that are meant for a specific purpose. Each of these parts have internal state variables that do not interact directly with the outside world. The interaction between components is made by exposing a visible interface that hides the internal complexity of the subassembly, and only displays linking parts. In programming this is called encapsulation. In a material world we call this a cover, or a panel. As an example consider an electricity panel. From the outside it only presents a groups of buttons, so the users of the panel will be able to control the system by interacting with them, and they will never get to play with the connections inside. Encapsulation in .NET is introduced by declaring either private or public flags in fields, properties or methods. Public declarations give your objects or modules an outside look, keeping the functional code hidden and protected.

Structure of the code

Lets start by saying that OpenVOGEL divides the code into two different branches: the calculation model, and the visual model. When you work with your project from the HMI (human machine interface), you are actually dealing with the visual model. When you launch a calculation, this model is internally converted into a calculation model, analysed and the converted back into another visual model for post processing. The reason why this has been done is simple: the calculation model does not need to know about how you represent your 3D model, and the visual model does not need to know nothing about how the results came out. Also, the calculation model has to deal with performance issues and does not need to know the details about how the geometry was generated. So the division has to do with "keeping on each side what is strictly necessary". Otherwise we would be bouncing all the time with irrelevant data, and when building a complex program, this can mean chaos, confusion and bugs.

In the next section we will describe the structure of the visual model contained in OpenVOGEL.DesignTools. The calculation model in OpenVOGEL.AeroTools is more complex because it deals with the aerodynamic and structural algorithms, so it will be subject of a dedicated section. If you can't way to see it, then navigate here.

Visual model

The visual model has two goals: gathering input data and showing the results to the user. So the visual model is naturally divided in two parts: the DesignModel and the ResultModel. Through the DesignModel you are confronted with a set of tools to declare how you want your model to be. When you run any simulation, the DesignModel is converted into a calculation object, and during the simulation results files are written to the hard disk. After the analysis, these files are read and translated into a ResultModel which confronts you with the results. The ResultModel contains the original settings and a collection of ResultFrames (depending on the simulation type).

The visual model lets you create an airplane by assembling together a number of four different kind of objects:

• LiftingSurface
• Fuselage
• JetEngine
• ImportedSurface

The three standard surface types (wings, fuselage and engines) assembled together to represent an airplane model.

The results are displayed in a ResultModel

As said before, these four classes inherit a common ancestor called Surface. The DesignModel stores all the surfaces in a List(Of Surface), and provides public methods to add every particular type into the heap.

Public Class DesignModel
   [...]
   Public Property Objects As New List(Of Surface)
   Public Sub AddLiftingSurface()
   Public Sub AddExtrudedBody()
   Public Sub AddJetEngine()
   [...]
End Class

Because of the surfaces are all treated equally, you could easily add your own types to the software. Basically you would just create a new Surface descendant type, and write a method on the DesignModel to load it on the heap.

Public Class MyNewSurfaceType 
   Inherits Surface 
   [...]
End Class

Public Class DesignModel
   [...]
   Public Sub AddMySurfaceType() 
      Dim NewSurface = New MyNewSurfaceType
      NewSurface.Name = String.Format("Surface - {0}", Objects.Count) 
      Objects.Add(NewSurface) 
   End Sub 
   [...]
End Class

Of course you would also need to make a user control if you want HMI support for it, but the above code snippet would make the basic trick. If only working in the console project, you would definitely omit the HMI and simply implement the object directly in the code.

Project root

Now you know how the models are structured, you are probably wondering where the DesignModel and the ResultModel actually reside. Well, OpenVOGEL stores them on the ProjectRoot static module, which is located in the DesignTools/DataStore directory. You can access this module throughout the whole application. ProjectRoot provides not only access to most of the program data, but it also contains the logic that manages the three different user interface modes: design, calculation setup and post-processing of results. It does this by exposing a set of public calls:

• Calls for synchronous calculation startup and loading of results
• Calls for input/output of OpenVOGEL (*.vog) files.

So much of what you do on the HMI, is actually handled here.

Public Module ProjectRoot
   Public Property Name As String = "New aircraft"
   Public Property FilePath As String = "" 
   Public Property SimulationSettings As New SimulationSettings 
   Public Property Model As DesignModel 
   Public Property Results As New ResultModel 
   Public Property VelocityPlane As New VelocityPlane 
   Public Property CalculationCore As Solver
   [...]
End Module

Components

Lets now take a closer look at the three standard components the software currently provides. You have probably already noticed that the greatest difference between them resides in the way the mesh is created, rather than how the are handled. They all expose a set of special parametric properties that are intended to generate a particular kind of mesh, and they do this by overriding the GenerateMesh method that they inherit from Surface.

Public MustInherit Class Surface 
   [...]
   Public Overridable Sub GenerateMesh()
   [...]
End Class

Public Class LiftingSurface 
   Inherits Surface 
   [...]
   Public Overrides Sub GenerateMesh()
      [...]
   End Sub
End Class

Lifting surfaces

Probably the most popular type is the LiftingSurface class. This surface type has been equipped with a special meshing algorithm that lets you model slender wings by declaring a row of adjacent macro panels represented by the class WingRegion. This class gathers all of the properties that are necessary to describe a single intermediate region of the wing, like the tip chord, the length and the sweepback angle. It also contains the parameters that define the local mesh as a grid: the number of span-wise and chord-wise panels.

Public Class WingRegion
   [...]
   Public Property SpanPanelsCount As Integer
   Public Property ChordNodesCount As Integer
   Public Property TipChord As Double
   Public Property Length As Double
   Public Property Sweepback As Double
   [...]
End Class

Public Class LiftingSurface  
   [...]
   Public Property WingRegions As New List(Of WingRegion)
   [...]
End Class

Fuselages

The fuselage type has a radically different meshing technique. Fuselages in OpenVOGEL are a sort of lofted surfaces defined from a set of longitudinal cross sections that are interpolated to locate the mesh nodes.

The real complexity of the fuselages, however, does not reside here but in the necessity of wing anchors. These features are needed in the panel method to provide continuity in the circulation and to avoid leakage. Before generating the mesh nodes, the meshing algorithm has to scan all connected wings and generate the anchor lines by projecting the wing root nodes on the primitive fuselage surface. Once the anchors are ready the surface is meshed in longitudinal chunks. If a chunk contains an anchor the mesh nodes are forced to keep the interface points in position.

The wing anchors are generated with the local longitudinal axis coincident to the global X axis. This limitation makes the anchoring algorithm easier, but it is not compatible with post meshing transformations (translation/rotation/translation). Therefore, fuselages that are meant to be anchored should not be rotated or translated.

Example of anchor lines joining lifting surfaces to the fuselage.

Public Class Fuselage
   Inherits Surface
   [...]
   Public Property CrossSections As List(Of CrossSection)
   Public Property AnchorLines As List(Of AnchorLine)
   Public Property MeshType As MeshTypes = MeshTypes.StructuredQuadrilaterals
   [...]
End Class
End Class

Nacelles

In short, nacelles are simply thin wall tubes. These surfaces offer an option to shed a closed wake from their trailing edge, so they can be useful to analyse jet engines or fan-ducts. This is specially interesting to have an idea of how the lift is transferred from the wing to the engine.

Imported surfaces

Imported surfaces are the only kind of surface that are not created in a parametric way. They are not created internally, but just read from a file that has been created manually or by a third party program. The file is only read when editing the surface and then stored internally for further use. The original mesh is kept in the project XML file when saving the project, so that when the model is reopened the original file is no longer required.

These surfaces are part of our effort to create interoperability with other software. At the moment, APAME files can be imported after some modifications in a text editor. There is still no standard for the input files, but the source code can always be used as reference.

3D representation of the models in OpenGL

The representation of the models in Tucan is done through OpenGL calls. Due to historical reasons, OpenVOGEL curretly only uses compatibility mode calls (i.e. OpenGL 1.4) instead of Core Profile. As explained at the beginning of this chapter, OpenGL is linked to the project through the SharpGL binding written in C# by Dave Kerr. That project does not only provide the necessary OpenGL context, but also a Winforms widget that displays the resulting pixel buffer in a winforms application and exposes the related drawing events.

Normally one would expect to find the drawing procedures directly inside the different clases. However, because the idea of the project was to provide independent libraries for different purposes, the OpenGL dependencies have been located in Tucan itself (the only project with HMI). Neither AeroTools nor DesignTools contain reference to SharpGL. The drawing procedures have been located inside a single unit and written as extensions to the parent classes. Extending a class by adding extra procedures is a quite useful and easy to implement .NET feature.

The rendering procedures can be found in OpenVOGEL.Tucan.Utility.ModelRendering.vb. There is a dispatching procedure that can be used for any surface, and a particular procedure for each surface kind.

Module ModelRendering
    ''' General extension that redispatches the rendering method to the correct surface:
    <Extension()>
    Public Sub Refresh3DModel(This As Surface,
                              ByRef gl As OpenGL,
                              Optional ByVal ForSelection As Boolean = False,
                              Optional ByVal ElementIndex As Integer = 0)
        If TypeOf This Is LiftingSurface Then
            Dim Surface As LiftingSurface = This
            Refresh3DModel(Surface, gl, ForSelection, ElementIndex)
        ElseIf TypeOf This Is Fuselage Then
            Dim Surface As Fuselage = This
            Refresh3DModel(Surface, gl, ForSelection, ElementIndex)
        ElseIf TypeOf This Is JetEngine Then
            Dim Surface As JetEngine = This
            Refresh3DModel(Surface, gl, ForSelection, ElementIndex)
        ElseIf TypeOf This Is ImportedSurface Then
            Dim Surface As ImportedSurface = This
            Refresh3DModel(Surface, gl, ForSelection, ElementIndex)
        ElseIf TypeOf This Is ResultContainer Then
            Dim Surface As Fuselage = This
            Refresh3DModel(Surface, gl, ForSelection, ElementIndex)
        End If
    End Sub
    ...
    ''' All the extensions for the particular classes:
    <Extension()>
    Public Sub Refresh3DModel(This As LiftingSurface,
                              ByRef gl As OpenGL,
                              Optional ByVal ForSelection As Boolean = False,
                              Optional ByVal ElementIndex As Integer = 0)
    ...
    End Sub
    ...
...
End Module

The actual OpenGL signals are contained in OpenVOGEL.Tucan.Utility.ModelInterface.vb. The signals include:

• Rebuilding the OpenGL lists (OpenGL 1.4) for each surface when the model has changed.
• Refreshing the model using the lists and the current camera setup when the 3D container requests a redrawing.
• Handling the selection of the components using the OpenGL hits when the user clicks on the 3D container.
• Representing the transit state (animation) in the 3D container.

The aerodynamic model in OpenVOGEL

The calculation core (CC) is the part of the program that deals with the calculation algorithms. As explained in the chapter about the source code, OpenVOGEL has been written in such a way that the calculation core is independent from the visual model. The only moment the calculation model meets the visual model is at calculation startup, when the later is converted into the former. The calculation core is also purely based in object orientation, and the panel method is very well suited for this. In fact, every panel type has been implemented by a class. Similarly, the finite elements method for the structural part of the aeroelastic analysis has been programmed in an object oriented style. We will start this chapter by taking a look at the main components of the aerodynamic model. After this we will examine the mathematics of the aerodynamic problem in more details.

The aerodynamic calculation core in classes

OpenVOGEL structures the panel method through the introduction of two main basic types: the Vortex class and the VotexRing interface. The Vortex class, representing a straight vortex segment, and the VortexRing interface, representing flat panels. This later is implemented by two classes: the VortexRing3, representing a triangle, and the VortexRing4, representing a quadrilateral. All these element types are based in the important concept of node, which is represented by the Node class. Nodes are actually the only material points that are traced in the fluid domain.

Public Class Node
   [...]
   Public Position As Vector3
   Public Displacement As Vector3
   Public Velocity As Vector3
   [...]
End Class

Public Class Vortex
   [...]
   Public Node1 As Node
   Public Node2 As Node
   [...]
End Class

Public Interface VortexRing
   [...]
   Property G As Double
   Property Node(ByVal Index As Integer) As Node
   [...]
End Interface

Public Class VortexRing3
   Implements VortexRing
   Private _Nodes(2) As Node 
   [...]
End Class

Public Class VortexRing4
   Implements VortexRing
   Private _Nodes(3) As Node 
   [...]
End Class

One step higher in the categorization of objects we find the lattices, which are represented by the general Lattice class. A lattice is a collection of VortexRing and/or Vortex elements that interconnect a cloud of Node objects.

Public Class Lattice
   [...]
   Public Nodes As New List(Of Node) 
   Public VortexRings As New List(Of VortexRing) 
   Public Vortices As New List(Of Vortex) 
   [...]
   Public Sub AddInducedVelocity(ByRef Velocity As Vector3, ByVal Point As Vector3, ByVal CutOff As Double)
   [...]
End Class

This means that a lattice can contain a mix of triangular and quadrilateral panels, plus a network of interconnected vortices. Lattices expose public methods to calculate the velocity induced by its components at a given point in space, and these methods are used to build the aerodynamic influence coefficients and to shed wakes.

There are in the code two extra classes derived from the Lattice class. These are the Wake and the BoundedLattice classes, which have very specific meanings. A Wake is a lattice that represents the free vortex sheets, while a BoundedLattice is a lattice that represents a solid surface (or more appropriate to say, a thin boundary layer).

Public Class Wake
   Inherits Lattice
   Public Property Primitive As New Primitive
   Public Property CuttingStep As Integer
   Public Property SupressInnerCircuation As Boolean
   [...]
End Class

Public Class BoundedLattice
   Inherits Lattice
   Public Wakes As New List(Of Wake)
   Public Sub PopulateWakeRings(Dt As Double, TimeStep As Integer)
   Public Sub PopulateWakeRingsAndVortices(Dt As Double, TimeStep As Integer)
   [...]
End Class

Bounded lattices contain a list of wakes, and wakes contain a reference to the primitive edge panels on the bounded lattices from which they are shed into the fluid domain. These primitive panels provide the value of circulation that will characterize the wake rings and or vortices during its complete lifetime. To declare a primitive you have to introduce the indices of the nodes at the shedding edge, and then the indices of the corresponding panels, both in consecutive order.

Public Class Primitive
   Public Nodes As New List(Of Integer)
   Public Rings As New List(Of Integer)
End Class

The definition of a primitive edge can clearly be seen on the next picture. It is also clear from this picture that the primitive edge is not bounded to the root and the tip of the wing, but it can be defined anywhere, including the leading edge.

Definition of a primitive shedding edge.

Storage of wake nodes.

Wake lifetime and special modeling features

A remarkable feature of the calculation core is that every wake can adopt its own lifetime. This allows the user to control the length of each wake, and also the boundary conditions. With a normal wake, for instance, a problem arises in the wing root exactly at the point the wing trailing edge is merged with the fuselage. Shedding a wake node from there is not a good idea because spurious velocity components are introduced in the proximity of the surfaces. Still we can impose the Kutta condition in this small portion of the trailing edge without any extra code by introducing a wake of 0 CuttingStep. OpenVOGEL does this automatically in the conversion module, and the kernel does not know anything about the real use of that wake.

A different issue that does require the kernel to be aware of is the suppression of the inner circulation. This problem arises again at the wing root, but in the direction of the free stream. The vortices that are adjacent to the fuselage are not so free to move in reality, since they are very close to the fuselage. Therefore, to avoid rolling the inner part of the wake, the first wake vortex line at the wing root in the direction of the flow must be suppressed. The kernel can do this for each wake by activating the SupressInnerCircuation property.

Building blocks: vortices, doublets and sources

Now we have introduced the main classes of the aerodynamic calculation core, we can go more into details. But before start talking about the main calculation algorithms, we probably need to understand what the vortices and vortex rings are capable of. As said before, VortexRing derived classes are found in two flavors, of three and four nodes respectively. But any of these two classes can also behave in two different ways: as flat constant-doublet or as constant-source distribution regions. For this purpose, the two classes implement methods that compute the doublet velocity, the doublet potential, the source velocity, and the source potential at any given point in space. The next set of functions encapsulated inside the VortexRing-derived classes constitute fundamental resources in the calculation core, as they allow determining the aerodynamic influence coefficients.

Public Interface VortexRing
   [...]

   ' Doublet velocity: adds the doublet velocity influence of the panel at the given Point to Vector
   Sub AddDoubletVelocityInfluence(Vector As Vector3, Point As Vector3, CutOff As Double, WithG As Boolean)

   ' Doublet potential: returns the doublet potential of the panel at the given Point
   Function GetDoubletPotentialInfluence(Point As Vector3, WithG As Boolean) As Double

   ' Source velocity: adds the source velocity influence of the panel at the given Point to Vector
   Sub AddSourceVelocityInfluence(Vector As Vector3, Point As Vector3, CutOff As Double, WithS As Boolean)

   ' Source potential: returns the source potential of the panel at the given Point
   Function GetSourcePotentialInfluence(Point As Vector3, WithS As Boolean) As Double

   [...]
End Interface

When a panel behaves as a constant doublet distribution, then the first two methods are used. The velocity influence is useful to calculate the induced velocity at any given point in space. If the argument WithG is set to true, then the local value of circulation is used, and the result is an absolute velocity. If WithG is set to false, then a unit circulation is used, and the result is the induced velocity per unit of circulation. Under the same behavior, the second function returns the value of the potential function at any given point in space, and the argument WithG plays the same role as before. When the panel behaves as a constant source distribution, the third function returns the value of the potential function at any point in space. The argument WithS is similar to WithG, indicating if a unit source/sink strength must be used.

You are now probably wondering in which situations we use either the doublet behavior or the source behavior of a panel, and when do we need to compute either the induced velocity or the local potential. Well, the answer to this is that it depends on the imposed boundary conditions. On slender surface panels we will always impose Neumann boundary conditions, while over thick surface panels we will always impose Dirichlet boundary conditions. So the whole boundary condition problem is naturally divided in two.

When the Neumann boundary condition is impose on a panel, the induced velocity must be computed at its middle control point located right at the surface. So we will scan all the surface and wake panels and use their doublet and source behavior to compute the local induced velocity. On the other hand, when the Dirichlet boundary condition is imposed on a panel, the potential must be computed at the inner control point located immediately under the surface, and for this we scan all the surface and wake panels to compute the velocity potential.

Local coordinates

Projection of a quad panel into a plane.

In OpenVOGEL the VortexRing derived classes implement the calculation of the potential functions in a local coordinates system. For this, a reference basis of orthogonal vectors $\mathbf {u} ,\mathbf {v} ,\mathbf {w}$ is created. For quadrilateral panels, $\mathbf {u}$ points in the direction of the first diagonal (from node 1 to node 3), $\mathbf {w}$ is normal to both diagonals, and $\mathbf {v}$ is normal to the other two vectors. For triangular panels, $\mathbf {u}$ points in the direction of the first edge (from node 1 to node 2), $\mathbf {w}$ is normal to the surface (which is always flat), and $\mathbf {v}$ is normal to the other two vectors. In all cases, vector $\mathbf {v}$ is always chosen so that the basis is dextrorotation. When a panels is declared as reversed, all vectors are just flipped.

In the following sections we will be using the notation $(u,v,w)$ to refer to coordinates in this local basis and respect to the panel control point. Most of the bibliography uses $(x,y,z)$ for this, but it is confusing to use the same notation as for the global coordinates system (referring to direction vectors $\mathbf {i} ,\mathbf {j} ,\mathbf {k}$). When we refer to velocity components, we will use $(V_{x},V_{y},V_{z})$ or $(V_{u},V_{v},V_{w})$. In the code, however, we are forced to use the letters x, y and z, since we are using in all cases the same Vector3 class.

Using this local basis, the position of the vertices is recomputed in two dimensions and cached for improved performance, so that they do not have to be recalculated continuously during execution (as it is done, for instance, in the FORTRAN examples of Katz & Plotkin). For the same reason, the length of the edges is also cached.

Note that the declaration of the local basis is not a requisite of the VortexRing interface. Actually, any implementation of the interface is free to choose its own projection method in its private part. In the current version of the CC the local basis method has also been implemented in the triangles regardless on the fact that they are always flat.

Collocation points

As explained before, the boundary condition problem requires the declaration of several collocation points. In OpenVOGEL, three collocation points are computed for each panel: the inner control point, the middle control point, and the outer control point. The inner one is only used to compute the inner potential of the thick bodies. The middle one is used as origin for the local coordinates (to complete the definition of the projected panel) and to compute the velocity in slender surfaces. Finally, the outer one is used to evaluate the velocity on thick bodies. The position of the inner and outer points is computed using the middle point and the normal direction times an epsilon scalar. The middle point is always computed using the average coordinates of the panel nodes.

Vortices

Although the whole calculation core is basically based on VortexRings (panels), there is a reason why it is sometimes useful to work with a more basic type: the Vortex. This situation is found when shedding wakes and when Neumann boundary conditions need to be imposed, that is to say, when there are slender surfaces in the model. Under this situation the induced velocity needs to be computed at the control points of the slender panels and at the nodal points of the wakes, and for this it is more efficient to use Vortices. In fact, when we use a VortexRing to compute the velocity somewhere, we loop around its three or four boundary segments. If we do this for two adjacent rings, then we will visit their common edge segment twice. But by implementing vortices we can avoid this and reduce the number of computations to almost the half. So the main reason why we use the Vortex element is because of computational efficiency.

The best way to take advantage of this is to equip the lattices with two parallel structures, a lattice of VortexRings and a lattice of Vortices, and to use one or the other in the most convenient situation. In OpenVOGEL the calculation core will always work with VortexRings and a parallel lattice of Vortices in the bounded lattices. On Wakes the situation is different: it will always work with a lattice of Vortices, and VortextRings will only be used if Dirichlet boundary conditions need to be applied somewhere (because they are needed for the computation of the potential). By doing this we have double gain: we reduce the size of the matrix problem by using VortexRings because there are always less VortexRings than Vortices, and we reduce the number of computations needed to build the right hand side and shed the wakes by using Vortices because they require a smaller number of operations.

When we work with a parallel lattice of vortices it is very efficient to give each Vortex a reference to the two or three adjacent VortexRings, along with information about their sense (that says how to interpret the sign of their circulation). With all this data, the intensity of each vortex can be evaluated by a simple sum. The process of generating the vortices of a lattice based on the VortexRings and searching and assigning the adjacent VortexRings to each Vortex is called PopulateVortices. This method is located inside the Lattice class (so it is common to all lattices). Note that the process does not need to scan all the lattices, but must be able to find at least 3 adjacent rings (because this occurs at the wing-body anchors).

For wakes the situation is again different. Since the circulation of wake rings and vortices remains constant during their whole life-cycle, their circulation is assigned directly at creation time during the wake shedding process.

Global adjacency survey

The continuity of circulation along the interface connecting two different surfaces is taken into account thanks to the gloval adjacency survery process.

Simulation of an extended flap by a second surface.

To calculate the jump of pressure across a slender VortexRing it is necessary to compute the local vorticity vector by vector-summing the side vortices. This operation requires knowing which VortexRings are adjacent to each side of each VortexRing. The process in charge of this global adjacency survey is called FindSurroundingRingsGlobally and is located in the source file AeroTools.CalculationModel.Solver_Calculations. This process must scan all the lattices, because at some interfaces (such as at the wing-body anchors), the rings might not share the same nodes. Two nodes are considered as being at the same point when their distance is less than a given tolerance.

The global adjacency survey process guarantees that the jump of pressure across panels that are adjacent is correctly calculated even when the panels do not share the same nodes. The only requisite for this is that the nodes connecting the shared edges are sufficiently close to each other (closer than the SurveyTolerance parameter declared in Settings). A practical application for this feature is the modeling of fowler flaps or control surfaces. If this process would not be included in the calculation core, the continuity in the circulation along the two surfaces would be broken, and each panel would see an opening at the shared edge. The picture here next represents a practical example. Note how the shedding edges have been declared, and how the lift is increased behind the flap hinge line.

In addition to all this, the global adjacency survey process is also necessary for computing the local velocity on thick bodies using the circulation gradient (this is explained later in the section about the Dirichlet boundary conditions).

The math problem

Introduction

Because OpenVOGEL deals with several types of panels and boundary conditions, the math and the implemented algorithms are a bit more complex than those required when only working with vortex rings.

The basic problem when dealing with panels is to find out the circulation that results in no normal flow across their surface. For slender panels this requirement is fulfilled by stating that the normal component of the fluid velocity at the control points must equal null. This is the so called Neumann boundary condition. Because the velocity associated to a vortex ring at a given point depends linearly on the circulation of that ring, a system of linear equations can be written, so that when solved, it will provide the value of the circulations that will cancel the normal velocity at every control point.

Now, when closed bodies such as fuselages need to be modeled, the previous way of proceeding produces sometimes an ill conditioned system that cannot be solved. The work around is to include sink/source panels and to impose the non-penetration condition by stating that total potential inside the body must equal a given constant value (typically chosen as zero). This is the so called Dirichlet boundary condition. Note that this is a very different kind of problem that requires calculating the potential associated to doublet and sink/source distributions. Additionally, in order to evaluate the surrounding flow velocity field, we will also have to calculate the velocity influence of these sink/source panels afterwards (the local surface velocity is approached differently, as explained later). So in the most general case, four basic algorithms are required:

• Calculation of the velocity of a vortex ring
• Calculation of the potential of a flat distribution of doublets
• Calculation of the velocity of a flat distribution of sink/sources
• Calculation of the potential of a flat distribution of sink/sources

So before we can formulate the complete boundary conditions and solve for the circulation of the panels, it is necessary to write down the mathematical expressions for potential and potential derivatives (velocity) at any point in space associated to a panel for both, a uniform sink/source distribution, and a uniform doublets distribution.

Note that we are under potential flow, so the next relationship holds:

$\nabla \phi (\mathbf {R} )=\mathbf {V} (\mathbf {R} )$

Doublet panels: velocity

Doublet panels represent a vortex loop of three or four straight vortex segments. The contribution of a panel to the velocity at any point in space is given by the celebrated Biot-Savart formula. The integrated equation yields:

$\mathbf {V} _{G}(\mathbf {R} )={\frac {G}{4\pi }}\sum _{k=1}^{n}{\frac {\mathbf {L} _{ji}\times \mathbf {r} _{i}}{||\mathbf {L} _{ji}\times \mathbf {r} _{i}||^{2}}}[\mathbf {L} _{ji}.(\mathbf {e} _{i}-\mathbf {e} _{j})]$

The summation extends to each side of the panel (i.e. $k=3$ for triangles, and $k=4$ for quadrilaterals). The sub-indices (i,j) refer to the first and second nodes of the side $k$, so they adopt the values $(1,2);(2,3);(3,1)$ for triangles and $(1,2);(2,3);(3,4);(4,1)$ for quadrilaterals. Vectors $\mathbf {L} _{ji}$ are the side segments given by $\mathbf {L} _{ji}=\mathbf {R} _{j}-\mathbf {R} _{i}$, where $\mathbf {R} _{i}$ and $\mathbf {R} _{j}$ are the position vectors of the two nodes on side $k$. Vector $\mathbf {r} _{i}$ is the relative position of the targeting point $\mathbf {R}$ respect to $\mathbf {R} _{i}$, that is to say, $\mathbf {r} _{i}=\mathbf {R} -\mathbf {R} _{i}$, and the unit vectors $\mathbf {e} _{i}$ and $\mathbf {e} _{j}$ point from each respective side node to the targeting point.

The VB.NET numerical implementation of the vortex segment velocity is as follows.

''' <summary>
''' Calculates BiotSavart vector at a given point. If WidthG is true vector is scaled by G.
''' </summary>
''' <remarks>
''' Calculation has been optimized by replacing object subs by local code.
''' Value types are used on internal calculations (other versions used reference type EVector3).
''' </remarks>
Public Function InducedVelocity(ByVal Point As Vector3,
                                Optional ByVal CutOff As Double = 0.0001,
                                Optional ByVal WithG As Boolean = True) As Vector3
    Dim Vector As New Vector3
    Dim D As Double
    Dim F As Double
    Dim Lx, Ly, Lz As Double
    Dim R1x, R1y, R1z, R2x, R2y, R2z As Double
    Dim vx, vy, vz As Double
    Dim dx, dy, dz As Double
    Dim NR1 As Double
    Dim NR2 As Double
    Lx = Node2.Position.X - Node1.Position.X
    Ly = Node2.Position.Y - Node1.Position.Y
    Lz = Node2.Position.Z - Node1.Position.Z
    R1x = Point.X - Node1.Position.X
    R1y = Point.Y - Node1.Position.Y
    R1z = Point.Z - Node1.Position.Z
    vx = Ly * R1z - Lz * R1y
    vy = Lz * R1x - Lx * R1z
    vz = Lx * R1y - Ly * R1x
    D = FourPi * (vx * vx + vy * vy + vz * vz)
    If D > CutOff Then
        ' Calculate the rest of the geometrical parameters:
        R2x = Point.X - Node2.Position.X
        R2y = Point.Y - Node2.Position.Y
        R2z = Point.Z - Node2.Position.Z
        NR1 = 1 / Math.Sqrt(R1x * R1x + R1y * R1y + R1z * R1z)
        NR2 = 1 / Math.Sqrt(R2x * R2x + R2y * R2y + R2z * R2z)
        dx = NR1 * R1x - NR2 * R2x
        dy = NR1 * R1y - NR2 * R2y
        dz = NR1 * R1z - NR2 * R2z
        F = (Lx * dx + Ly * dy + Lz * dz) / D
        If WithG Then
            F *= G
        End If
        Vector.X += F * vx
        Vector.Y += F * vy
        Vector.Z += F * vz
    Else
        Vector.X += 0
        Vector.Y += 0
        Vector.Z += 0
    End If
    Return Vector
End Function

Doublet panels: potential

Without going into details, the potential of a flat panel due to a uniform doublet distribution can by obtained from :

$\phi _{G}(\mathbf {R} )={\frac {G}{4\pi }}\sum _{k=1}^{n}[tan^{-1}({\frac {m_{ij}e_{i}-h_{i}}{zr_{i}}})-tan^{-1}({\frac {m_{ij}e_{j}-h_{j}}{zr_{j}}})]$

This formulation is written in local coordinates and the summation extends to each side of the panel (i.e. $k=3$ for triangles, and $k=4$ for quadrilaterals). The sub-indices (i,j) refer to the first and second nodes of the side $k$, so they adopt the values $(1,2);(2,3);(3,1)$ for triangles and $(1,2);(2,3);(3,4);(4,1)$ for quadrilaterals. Note that the above equation is further simplified in the numeric algorithms to avoid calculating two arc-tangents and gain performance.

The formulation is completed with the next definitions:

$m_{ij}={\frac {v_{j}-v_{i}}{u_{j}-u_{i}}}$

$r_{i}={\sqrt {(u-u_{i})^{2}+(v-v_{i})^{2}+w^{2}}}$

$e_{i}=(u-u_{i})^{2}+w^{2}$

The VB.NET numerical implementation of the above equation is as follows.

''' <summary>
''' Returns the influence of the doublet distribution in the velocity potential.
''' </summary>
Public Function GetDoubletPotentialInfluence(ByVal Point As Vector3,
                                             Optional ByVal WithG As Boolean = True) As Double Implements VortexRing.GetDoubletPotentialInfluence
    ' Convert the point to local coordinates (center on the control point and using the local basis)
    Dim dx = Point.X - _MidleControlPoint.X
    Dim dy = Point.Y - _MidleControlPoint.Y
    Dim dz = Point.Z - _MidleControlPoint.Z
    Dim p As New Vector3(dx * _Basis.U.X + dy * _Basis.U.Y + dz * _Basis.U.Z,
                         dx * _Basis.V.X + dy * _Basis.V.Y + dz * _Basis.V.Z,
                         dx * _Basis.W.X + dy * _Basis.W.Y + dz * _Basis.W.Z)
    Dim pzsq = p.Z * p.Z
    ' Distances:
    Dim r0px = p.X - _LocalNodes(0).X
    Dim r0py = p.Y - _LocalNodes(0).Y
    Dim r0p As Double = Math.Sqrt(r0px * r0px + r0py * r0py + pzsq)
    Dim r1px = p.X - _LocalNodes(1).X
    Dim r1py = p.Y - _LocalNodes(1).Y
    Dim r1p As Double = Math.Sqrt(r1px * r1px + r1py * r1py + pzsq)
    Dim r2px = p.X - _LocalNodes(2).X
    Dim r2py = p.Y - _LocalNodes(2).Y
    Dim r2p As Double = Math.Sqrt(r2px * r2px + r2py * r2py + pzsq)
    ' Use center point as reference to compute the altitude:
    Dim z As Double = p.Z
    ' Projected segments:
    Dim d01x As Double = _LocalEdges(0).X
    Dim d01y As Double = _LocalEdges(0).Y
    Dim d12x As Double = _LocalEdges(1).X
    Dim d12y As Double = _LocalEdges(1).Y
    Dim d20x As Double = _LocalEdges(2).X
    Dim d20y As Double = _LocalEdges(2).Y
    ' Entities for evaluation of arctangents:
    Dim z2 As Double = z * z
    Dim e0 As Double = r0px * r0px + z2
    Dim e1 As Double = r1px * r1px + z2
    Dim e2 As Double = r2px * r2px + z2
    Dim h0 As Double = r0px * r0py
    Dim h1 As Double = r1px * r1py
    Dim h2 As Double = r2px * r2py
    Dim f0 As Double = d01y * e0 - d01x * h0
    Dim g0 As Double = d01y * e1 - d01x * h1
    Dim tn01 As Double = Math.Atan2(z * d01x * (f0 * r1p - g0 * r0p), z2 * d01x * d01x * r0p * r1p + f0 * g0)
    Dim f1 As Double = d12y * e1 - d12x * h1
    Dim g1 As Double = d12y * e2 - d12x * h2
    Dim tn12 As Double = Math.Atan2(z * d12x * (f1 * r2p - g1 * r1p), z2 * d12x * d12x * r1p * r2p + f1 * g1)
    Dim f2 As Double = d20y * e2 - d20x * h2
    Dim g2 As Double = d20y * e0 - d20x * h0
    Dim tn20 As Double = Math.Atan2(z * d20x * (f2 * r0p - g2 * r2p), z2 * d20x * d20x * r2p * r0p + f2 * g2)
    Dim Potential As Double = (tn01 + tn12 + tn20) / FourPi
    If WithG Then Potential *= G
    Return Potential
End Function

Sink/source panels: potential

The calculation of the velocity potential associated to constant sink/source panels and vortex rings (constant doublets) at a given point in space is quite complex. OpenVOGEL reduces the complexity level by projecting the quadrilateral panels on its normal direction (according to the local basis), so that they can be represented by a flat panel (triangular panels do not need to be projected since they are always flat, although they are treated similarly). The calculation of the velocity potential under such situation can be found in this reference .

$\phi _{S}(\mathbf {R} )=-{\frac {S}{4\pi }}[\sum _{k=1}^{n}{\frac {(u-u_{i})(v_{j}-v_{i})-(v-v_{i})(u_{j}-u_{i})}{d_{ij}}}ln({\frac {r_{i}+r_{j}+d_{ij}}{r_{i}+r_{j}-d_{ij}}})-|w|\sum _{k=1}^{n}tan^{-1}({\frac {m_{ij}e_{i}-h_{i}}{wr_{i}}})-tan^{-1}({\frac {m_{ij}e_{j}-h_{j}}{wr_{j}}})]$

Of course, this method has some associated leakage, which depends on how well panels are approach by their flat counterpart. It is therefore important to provide a mesh where the quad panels are not too twisted.

When working with sink/source panels, the intensity associated to them is not an unknown. In fact, when the internal potential is targeted to zero, the sink/source intensity is just set to normal component of the free air-stream velocity:

$S_{i}=-\mathbf {V} _{\infty }\cdot \mathbf {n} _{i}$

The VB.NET numerical implementation of the potential for a triangular panel is as follows.

''' <summary>
''' Returns the influence of the source distribution in the velocity potential.
''' </summary>
Public Function GetSourcePotentialInfluence(ByVal Point As Vector3, Optional ByVal WithS As Boolean = True) As Double Implements VortexRing.GetSourcePotentialInfluence
    ' Convert the point to local coordinates (center on the control point and using the local basis)
    Dim dx = Point.X - _MidleControlPoint.X
    Dim dy = Point.Y - _MidleControlPoint.Y
    Dim dz = Point.Z - _MidleControlPoint.Z
    Dim p As New Vector3(dx * _Basis.U.X + dy * _Basis.U.Y + dz * _Basis.U.Z,
                         dx * _Basis.V.X + dy * _Basis.V.Y + dz * _Basis.V.Z,
                         dx * _Basis.W.X + dy * _Basis.W.Y + dz * _Basis.W.Z)
    Dim pzsq = p.Z * p.Z
    ' Distances:
    Dim r0px = p.X - _LocalNodes(0).X
    Dim r0py = p.Y - _LocalNodes(0).Y
    Dim r0p As Double = Math.Sqrt(r0px * r0px + r0py * r0py + pzsq)
    Dim r1px = p.X - _LocalNodes(1).X
    Dim r1py = p.Y - _LocalNodes(1).Y
    Dim r1p As Double = Math.Sqrt(r1px * r1px + r1py * r1py + pzsq)
    Dim r2px = p.X - _LocalNodes(2).X
    Dim r2py = p.Y - _LocalNodes(2).Y
    Dim r2p As Double = Math.Sqrt(r2px * r2px + r2py * r2py + pzsq)
    ' Use center point as reference to compute the altitude:
    Dim z As Double = Math.Abs(p.Z)
    ' Projected segments:
    Dim d01x As Double = _LocalEdges(0).X
    Dim d01y As Double = _LocalEdges(0).Y
    Dim d12x As Double = _LocalEdges(1).X
    Dim d12y As Double = _LocalEdges(1).Y
    Dim d20x As Double = _LocalEdges(2).X
    Dim d20y As Double = _LocalEdges(2).Y
    ' Segments length:
    Dim d01 As Double = _LocalSides(0)
    Dim d12 As Double = _LocalSides(1)
    Dim d20 As Double = _LocalSides(2)
    ' Logarithms:
    Dim ln01 As Double = (r0px * d01y - r0py * d01x) / d01 * Math.Log((r0p + r1p + d01) / (r0p + r1p - d01))
    Dim ln12 As Double = (r1px * d12y - r1py * d12x) / d12 * Math.Log((r1p + r2p + d12) / (r1p + r2p - d12))
    Dim ln20 As Double = (r2px * d20y - r2py * d20x) / d20 * Math.Log((r2p + r0p + d20) / (r2p + r0p - d20))
    ' Entities for evaluation of arctangents:
    Dim z2 As Double = z * z
    Dim e0 As Double = r0px * r0px + z2
    Dim e1 As Double = r1px * r1px + z2
    Dim e2 As Double = r2px * r2px + z2
    Dim h0 As Double = r0px * r0py
    Dim h1 As Double = r1px * r1py
    Dim h2 As Double = r2px * r2py
    Dim f0 As Double = d01y * e0 - d01x * h0
    Dim g0 As Double = d01y * e1 - d01x * h1
    Dim tn01 As Double = Math.Atan2(z * d01x * (f0 * r1p - g0 * r0p), z2 * d01x * d01x * r0p * r1p + f0 * g0)
    Dim f1 As Double = d12y * e1 - d12x * h1
    Dim g1 As Double = d12y * e2 - d12x * h2
    Dim tn12 As Double = Math.Atan2(z * d12x * (f1 * r2p - g1 * r1p), z2 * d12x * d12x * r1p * r2p + f1 * g1)
    Dim f2 As Double = d20y * e2 - d20x * h2
    Dim g2 As Double = d20y * e0 - d20x * h0
    Dim tn20 As Double = Math.Atan2(z * d20x * (f2 * r0p - g2 * r2p), z2 * d20x * d20x * r2p * r0p + f2 * g2)
    Dim Potential As Double = -(ln01 + ln12 + ln20 - z * (tn01 + tn12 + tn20)) / FourPi
    If WithS Then Potential *= S
    Return Potential
End Function

Sink/source panels: velocity

The velocity associated to a sink/source panel can be found in reference .

$V_{u}={\frac {S}{4\pi }}[\sum _{k=1}^{n}{\frac {v_{j}-v_{i}}{d_{ij}}}ln({\frac {r_{i}+r_{j}-d_{ij}}{r_{i}+r_{j}+d_{ij}}})]$

$V_{v}={\frac {S}{4\pi }}[\sum _{k=1}^{n}{\frac {u_{j}-u_{i}}{d_{ij}}}ln({\frac {r_{i}+r_{j}-d_{ij}}{r_{i}+r_{j}+d_{ij}}})]$

$V_{w}={\frac {S}{4\pi }}[\sum _{k=1}^{n}tan^{-1}({\frac {m_{ij}e_{i}-h_{i}}{zr_{i}}})-tan^{-1}({\frac {m_{ij}e_{j}-h_{j}}{zr_{j}}})]$

The velocity in the global coordinates system would be built by summing the three orthogonal projections:

$\mathbf {V} _{S}(\mathbf {R} )=V_{u}\mathbf {u} +V_{v}\mathbf {v} +V_{w}\mathbf {w}$

The VB.NET numerical implementation of the source velocity influence is as follows.

''' <summary>
''' Adds the influence of the source distribution in the velocity.
''' </summary>
''' <remarks></remarks>
Sub AddSourceVelocityInfluence(ByRef Vector As Vector3,
                               ByVal Point As Vector3,
                               Optional ByVal WithS As Boolean = True) Implements VortexRing.AddSourceVelocityInfluence
    ' Convert the point to local coordinates (center on the control point and using the local basis)
    Dim dx = Point.X - _MidleControlPoint.X
    Dim dy = Point.Y - _MidleControlPoint.Y
    Dim dz = Point.Z - _MidleControlPoint.Z
    Dim p As New Vector3(dx * _Basis.U.X + dy * _Basis.U.Y + dz * _Basis.U.Z,
                         dx * _Basis.V.X + dy * _Basis.V.Y + dz * _Basis.V.Z,
                         dx * _Basis.W.X + dy * _Basis.W.Y + dz * _Basis.W.Z)
    Dim pzsq = p.Z * p.Z
    ' Distances:
    Dim r0px = p.X - _LocalNodes(0).X
    Dim r0py = p.Y - _LocalNodes(0).Y
    Dim r0p As Double = Math.Sqrt(r0px * r0px + r0py * r0py + pzsq)
    Dim r1px = p.X - _LocalNodes(1).X
    Dim r1py = p.Y - _LocalNodes(1).Y
    Dim r1p As Double = Math.Sqrt(r1px * r1px + r1py * r1py + pzsq)
    Dim r2px = p.X - _LocalNodes(2).X
    Dim r2py = p.Y - _LocalNodes(2).Y
    Dim r2p As Double = Math.Sqrt(r2px * r2px + r2py * r2py + pzsq)
    ' Projected segments:
    Dim d01x As Double = _LocalEdges(0).X
    Dim d01y As Double = _LocalEdges(0).Y
    Dim d12x As Double = _LocalEdges(1).X
    Dim d12y As Double = _LocalEdges(1).Y
    Dim d20x As Double = _LocalEdges(2).X
    Dim d20y As Double = _LocalEdges(2).Y
    ' Segments length:
    Dim d01 As Double = _LocalSides(0)
    Dim d12 As Double = _LocalSides(1)
    Dim d20 As Double = _LocalSides(2)
    ' Use center point as reference to compute the altitude:
    Dim z As Double = Math.Abs(p.Z)
    ' Entities for evaluation of arctangents:
    Dim z2 As Double = z * z
    Dim e0 As Double = r0px * r0px + z2
    Dim e1 As Double = r1px * r1px + z2
    Dim e2 As Double = r2px * r2px + z2
    Dim h0 As Double = r0px * r0py
    Dim h1 As Double = r1px * r1py
    Dim h2 As Double = r2px * r2py
    ' Normal component:
    ' These are the Katz-Plotkin fotran formulas, which are based in only one Atan2 instead of two
    Dim f0 As Double = d01y * e0 - d01x * h0
    Dim g0 As Double = d01y * e1 - d01x * h1
    Dim tn01 As Double = Math.Atan2(z * d01x * (f0 * r1p - g0 * r0p), z2 * d01x * d01x * r0p * r1p + f0 * g0)
    Dim f1 As Double = d12y * e1 - d12x * h1
    Dim g1 As Double = d12y * e2 - d12x * h2
    Dim tn12 As Double = Math.Atan2(z * d12x * (f1 * r2p - g1 * r1p), z2 * d12x * d12x * r1p * r2p + f1 * g1)
    Dim f2 As Double = d20y * e2 - d20x * h2
    Dim g2 As Double = d20y * e0 - d20x * h0
    Dim tn20 As Double = Math.Atan2(z * d20x * (f2 * r0p - g2 * r2p), z2 * d20x * d20x * r2p * r0p + f2 * g2)
    Dim Vw As Double = Math.Sign(p.Z) * (tn01 + tn12 + tn20)
    ' Tangent components
    Dim ln01 As Double = Math.Log((r0p + r1p - d01) / (r0p + r1p + d01))
    Dim ln12 As Double = Math.Log((r1p + r2p - d12) / (r1p + r2p + d12))
    Dim ln20 As Double = Math.Log((r2p + r0p - d20) / (r2p + r0p + d20))
    ' Planar velocity componets:
    Dim Vu As Double = d01y / d01 * ln01 + d12y / d12 * ln12 + d20y / d20 * ln20
    Dim Vv As Double = d01x / d01 * ln01 + d12x / d12 * ln12 + d20x / d20 * ln20
    ' Recompose vector in global coordinates
    Dim Factor As Double = 1.0# / FourPi
    If WithS Then
        Factor *= S
    End If
    Vu *= Factor
    Vv *= Factor
    Vw *= Factor
    Vector.X += _Basis.U.X * Vu - _Basis.V.X * Vv + _Basis.W.X * Vw
    Vector.Y += _Basis.U.Y * Vu - _Basis.V.Y * Vv + _Basis.W.Y * Vw
    Vector.Z += _Basis.U.Z * Vu - _Basis.V.Z * Vv + _Basis.W.Z * Vw
End Sub

Note that at the end of the code the local velocity components (in local coordinates) are projected back to the global coordinates system using the local panel orthonormal basis.

Dirichlet boundary conditions

The so called Dirichlet boundary conditions are imposed on the surface of thick closed bodies, requiring the potential immediately under the surface to be constant. To write this in mathematical terms, we split the potential associated to bounded doublets and sources as follows:

$\phi _{G}(\mathbf {R} _{i})+\phi _{S}(\mathbf {R} _{i})+\phi _{W}(\mathbf {R} _{i})=constant=0$

In this expression, $\phi _{G}$ is the potential associated to the bounded doublet panels, $\phi _{S}$ is the potential associated to the bounded sink/source panels and $\phi _{W}$ is the potential associated to the wakes (which by definition, only contain doublets). The preference for making the constant equal zero is related to the previous assumption that sink/source intensities are determined by the normal component of the free stream velocity.

Since the sink/source intensity of each bounded panel is known in advance, for rigid panels (i.e. panels attached to nodes that do not displace) the unit potential at each bounded inner control point can be calculated once and stored in a source/sink potential matrix $\mathbf {M} _{S}$. This matrix can then be reused to compute the total sink/source potential when necessary as follows:

$\mathbf {\phi } _{S}=\mathbf {M} _{S}\mathbf {S}$

Note that this technique is maybe not so suitable for a lattice that is subject to deformations. However, in OpenVOGEL this problem is not encountered because the aeroelastic problem is reserved for slender panels only. So whenever we use thick bodies, the shape of the corresponding panels always remain intact during the whole calculation.

The unknowns of the Dirichlet problem are the doublet intensities of each panel surrounding the body, represented here by vector $\mathbf {G}$. Similarly to the sources, we can write a matrix of unit-intensity potential $\mathbf {M} _{G}$ so that:

$\mathbf {\phi } _{G}=\mathbf {M} _{G}\mathbf {G}$

Finally, the potential associated to the wakes needs to be added. Because wakes are constantly changing their shape, in place of building a matrix, we rather include their potential as a summation of individual contributions:

$\mathbf {\phi } _{Wi}=\sum _{j=1}^{N}\phi _{Wj}(\mathbf {R} _{i})$

Because wakes are shed from lifting surfaces, their potential is only associated to constant doublet panels of known intensity. When wake panels are far enough from the body, the general practice is to approach their contributions to the potential by point doublets, which have much simpler mathematical expressions. This solution is known as the far field potential. However, the current version of the CC does not include this feature yet. The complete equations are always used.

The Dirichlet problem is then represented by the next system or linear equations

$\mathbf {M} _{G}\mathbf {G} =-\mathbf {M} _{S}\mathbf {S} -\mathbf {\phi } _{W}$

which needs to be assembled and solved together along with the Neumann problem (see next section) in a global matrix problem.

Local surface velocity in thick closed bodies

Once the circulations $\mathbf {G}$ are found, the local velocity is not calculated by summing the velocity influences. A practical implementation demonstrates that that kind of approach does not yield correct results, because it is unable to resolve the local tangent velocity due to the surface circulation. Instead, the local velocity is approached through the local two dimensional circulation gradient. In general we would have:

$\mathbf {V} =(\mathbf {V} _{\infty }.\mathbf {u} -dG/du)\mathbf {u} +(\mathbf {V} _{\infty }.\mathbf {v} -dG/dv)\mathbf {v}$

The problem here is how to determine the local derivatives $dG/du$ and $dG/dv$. Of course there is no analytic formula for the function $G(u,v)$, since the circulation on each panel is assumed constant. So we need to use the information of the adjacent panels.

Tucan solves this in a relatively simple way by using least squares to fit a differential function using the adjacent panels circulation. The method is inspired in the central difference method (and the mean-value theorem), but the nice thing about this approach is that it works for any kind of panel without the slightest modification, as long as there are at least two or more none co-linear adjacent panels (which for a closed body should always be the case). The algorithm is:

1. Using the three or four adjacent panels and the panel itself, build matrix $\mathbf {M}$ as follows:

$\mathbf {M} =\left[{\begin{array}{ccc}\sum u_{i}.u_{i}&\sum u_{i}.v_{i}&\sum u_{i}\\\sum u_{i}.v_{i}&\sum v_{i}.v_{i}&\sum v_{i}\\\sum u_{i}&\sum v_{i}&n\end{array}}\right]$

Where $n$ is the total number of sample points (normally 4 or 5 for triangular and quadrilateral panels respectively).

2. Build right hand side vector $\mathbf {B}$ as follows:

$\mathbf {B} =\left[{\begin{array}{ccc}\sum u_{i}.G_{i}\\\sum v_{i}.G_{i}\\\sum G_{i}\end{array}}\right]$

3. Solve the 3x3 system of linear equations $\mathbf {M} \mathbf {A} =\mathbf {B}$.

4. Finally we will have:

$dG/du=A_{1}$

$dG/dv=A_{2}$

The third component of the solution vector is a mean circulation value, and can be discarded.

The next VB.NET code shows how this algorithm has been implemented:

Public Sub CalculateLocalVelocity(StreamVelocity As Vector3) Implements VortexRing.CalculateLocalVelocity
    VelocityT.SetToCero()
    If Not _IsSlender Then
        Dim M As New Matrix(3)
        Dim B As New Vector(3)
        ' Local circulation at (0,0)
        M(2, 2) += 1
        B(2) += G
        ' Add circulation of adjacent panels
        For i = _SurroundingRings.GetLowerBound(0) To _SurroundingRings.GetUpperBound(0)
            ' Do not include the panels behind a shared edge
            If _SurroundingRings(i, 0) IsNot Nothing AndAlso _SurroundingRings(i, 1) Is Nothing Then
                Dim Delta As New Vector3
                Delta.X = _SurroundingRings(i, 0).ControlPoint.X - ControlPoint.X
                Delta.Y = _SurroundingRings(i, 0).ControlPoint.Y - ControlPoint.Y
                Delta.Z = _SurroundingRings(i, 0).ControlPoint.Z - ControlPoint.Z
                Dim OtherU As Double = Delta.InnerProduct(_Basis.U)
                Dim OtherV As Double = Delta.InnerProduct(_Basis.V)
                Dim OtherG = _SurroundingRings(i, 0).G
                M(0, 0) += OtherU * OtherU
                M(0, 1) += OtherU * OtherV
                M(0, 2) += OtherU
                M(1, 0) += OtherU * OtherV
                M(1, 1) += OtherV * OtherV
                M(1, 2) += OtherV
                M(2, 0) += OtherU
                M(2, 1) += OtherV
                M(2, 2) += 1
                B(0) += OtherU * OtherG
                B(1) += OtherV * OtherG
                B(2) += OtherG
            End If
        Next
        ' Calculate the circulation derivatives on each tangent directions
        ' Vector A will contain the circulation slopes and the local mean circulation
        Dim Equations As New LinearEquations
        Dim A As Vector
        Try
            A = Equations.Solve(M, B)
        Catch ex As Exception
            A = New Vector(3)
        End Try
        ' Recompose velocity in global coordinates
        Dim StreamU As Double = _Basis.U.InnerProduct(StreamVelocity)
        Dim StreamV As Double = _Basis.V.InnerProduct(StreamVelocity)
        VelocityT.Add(_Basis.U, -A(0) + StreamU)
        VelocityT.Add(_Basis.V, -A(1) + StreamV)
    End If
End Sub

Note that when a panel has more than two adjacent panels in one side, that side is disregarded because it is assumed that one of the panels is blocking the other one (which is the case of a slender panel connected to the body as an anchor). An additional requirement for this method is that all adjacent panels in the model must be found (using the global adjacency survey process described before).

It might sound strange that the above numerical recipe actually produces the required velocity. However, as explained before, the trick is in using the circulation of the adjacent panels to best fit a two dimensional differential (i.e. linear) function centered on the target panel. The limitation of this algorithm might be in the curvature of the surface. Note that the gradient is approached using the projection of each relative position vectors on each local direction, thus it takes a shortcut to the adjacent control points instead of using real surface coordinates.

Validation of the algorithm on a sphere (for which the analytical solution is known) demonstrates however that the resulting local velocity and $C_{P}$ are predicted with very good accuracy (see picture here next).

Distribution of the pressure coefficient over the surface of an sphere (theory .vs. OpenVOGEL Tucan). The model in Tucan consisted of 2400 quadrilateral panels.

Pressure coefficient in thick bodies

Once the local surface velocity is determined, the pressure coefficient can be calculated using Bernoulli's equation. For the steady case, we have:

$C_{p}=1-({\frac {V}{V_{\infty }}})^{2}$

Cp = 1 - ((VelocityT.X * VelocityT.X + VelocityT.Y * VelocityT.Y + VelocityT.Z * VelocityT.Z)) / Vsqr

Neumann boundary conditions

In slender surfaces we apply the so called Neumann boundary conditions, which are expressed by stating that the normal velocity component across the surface must be null (i.e. no airflow across the surface). We can the total velocity at a control point in different components, and project it on the panel normal vector to obtain a net cross flow:

$Q_{i}=[\mathbf {V} _{\infty }+\mathbf {V} _{G}(\mathbf {R} _{i})+\mathbf {V} _{S}(\mathbf {R} _{i})+\mathbf {V} _{W}(\mathbf {R} _{i})]\cdot \mathbf {n} _{i}=0$

This equation must be repeated for each slender panel. We will name $N_{G}$ to the total number of slender panels, meaning that there will be a system of $N_{G}\times N_{G}$ equations. In a similar way as it has been done for the Dirichlet boundary conditions, here we can build an influence cross-flow matrix for the doublets, another one for the sources, and separate the wake cross-flow in a summation apart. As explained in the previous chapters, the velocity associated to the wake $\mathbf {V} _{W}(\mathbf {R} _{i})$ can better be calculated using the $N_{WV}$ wake vortices. We can write a wake related cross flow vector as follows:

$\mathbf {Q} _{W}=[\mathbf {Q} _{1},...,\mathbf {Q} _{i},...,\mathbf {Q} _{N_{G}}]$

with each component as follows:

$Q_{Wi}=\sum _{k=1}^{N_{WV}}\mathbf {V} _{Vk}(\mathbf {R} _{i})\cdot \mathbf {n} _{i}$

By $\mathbf {V} _{Vk}$ we mean here the influence velocity of the wake vortex $k$ obtained by the Biot-Savart equation using the vortex intensity that has been inherited from the primitive shedding edge.

The free stream cross-flow components can be written as follows:

$\mathbf {Q} _{\infty }=[\mathbf {V} _{\infty }\cdot \mathbf {n} _{1},...,\mathbf {V} _{\infty }\cdot \mathbf {n} _{i},...,\mathbf {V} _{\infty }\cdot \mathbf {n} _{N_{G}}]$

The sources cross-flow components can be written as follows:

$\mathbf {Q} _{S}=\mathbf {M} _{S}\mathbf {S}$

The matrix $\mathbf {M} _{S}$ contains the influence cross flow of each thick surface panel on each slender surface panel due to the thick surface panels sink/source strength. Therefore, the size of this matrix is $N_{G}\times N_{S}$ and each component can be written as follows:

$M_{Sij}=\mathbf {V} _{Sj}(\mathbf {R} _{i})\cdot \mathbf {n} _{i}$

By $\mathbf {V} _{Sj}$ we mean here the influence velocity associated to the flat sink-source distribution over panel $j$, which can be obtained from the formulas that were written in the previous paragraphs.

The doublets cross-flow components can be written as follows:

$\mathbf {Q} _{G}=\mathbf {M} _{G}\mathbf {G}$

The matrix $\mathbf {M} _{G}$ contains the influence cross flow of the slender surface panels due to the panels dipole strength (circulation). Therefore, the size of this matrix is $N_{G}\times N_{G}$ and each component can be written as follows:

$M_{Gij}=\mathbf {V} _{Gj}(\mathbf {R} _{i})\cdot \mathbf {n} _{i}$

By $\mathbf {V} _{Gj}$ we mean here the influence velocity associated to the flat dipole distribution over panel $j$ (vortex rings), which can be obtained from the formulas that were written in the previous paragraphs.

The unknowns of the problem are the dipole intensities (circulations) associated to the slender panels (vector $\mathbf {G}$), which can be found after solving the next system of linear equations:

$\mathbf {M} _{G}\mathbf {G} =-\mathbf {Q} _{\infty }-\mathbf {Q} _{S}-\mathbf {Q} _{W}$

This system needs to be merged with the Dirichlet problem in a global matrix problem.

Mixed boundary conditions

In OpenVOGEL it is allowed to model thick and slender surfaces simultaneously. To accomplish this, the previous sets of equations are merged to form a single system of linear equations.

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:

1. The air-loads from certain rings of the vortex lattice are sent to certain nodes on the structure.
2. 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.
3. 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.

A mechanic link

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.

A kinematic link

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

1. Calculate RHS
2. Build matrix
3. Calculate circulation
4. Calculate airloads
5. Transfer loads to structure
6. Predict modal motion with new load
7. Transfer motion to lattice
8. Calculate velocity on wake nodal points
9. Update wakes

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

1. If t > 1
   1. Calculate velocity on wake nodal points
   2. Update wakes
2. 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)

1. 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.

Free flight simulation in OpenVOGEL

One of the most recent features developed in the OpenVOGEL suit is the free flight simulation module. This module combines the unsteady aeroelastic solver (which provides air loads for a given flow field) with the numerical integration of the rigid body equations of motion (in 6 degrees of freedom). It resembles in some aspects the aeroelastic module, but here the flow field is directly manipulated from the outcome of the equations (in the aeroelastic module, the variation of the flow field is a consequence of the repositioning if the shedding edges).

Coupling algorithm

The free flight simulation consists of an algorithm that couples the air loads and the motion. The air-loads are first computed using the instantaneous flow field at the beginning of each time step using the previous state. Based on these air loads, the motion is predicted (through numerical integration) for the next time step and the result is translated into a new flow field. The air loads are then recomputed in the new state and the motion is corrected in a cycle that is repeated until convergence is reached. Once the convergence criteria are met, the time is advanced by one step.

Equations of motion

If we attempt to write the equations of motion of a rigid body in an inertial reference frame, the inertial properties of the body become time dependent. So to reduce the complexity of the equations, the motion can better be described from the point of view of the body itself (a non inertial reference frame), where the inertia tensor is an invariant. Additionally, to avoid having the kinetic moment being coupled to the momentum, it is convenient to take the reference point for the moments as the center of gravity. Finally, to avoid a complex coupling of the angular velocities in the angular momentum equations, it is much simpler to align the axes with the main axes of inertia.

With all this said, the equations of motion for the rigid aircraft model including the external forces and gravity field can be written as follows :

${\frac {^{R}d}{dt}}{\mathbf {G}}={\mathbf {F}}+m{\mathbf {g}}$

${\frac {^{R}d}{dt}}{\mathbf {H}}={\mathbf {M}}$

Note that we have used Tenembaum's notation, in which the superscript R at the left of the differential operator means that the derivatives are taken with respect to the inertial reference frame (for the sake of simplicity, the earth in our case).

To convert the derivatives to the non inertial reference frame A (the moving aircraft), we use the reference frame transformation formula:

${\frac {^{R}d}{dt}}{\mathbf {a}}={\frac {^{A}d}{dt}}{\mathbf {a}}+{\mathbf {\omega }}\times {\mathbf {a}}$

where ${\mathbf {a}}$ is a generic vector.

The equations turn now to be:

${\frac {^{A}d}{dt}}{\mathbf {v}}+{\mathbf {\omega }}\times {\mathbf {v}}={\frac {1}{m}}{\mathbf {F}}+{\mathbf {g}}$

$II.{\frac {^{A}d}{dt}}{\mathbf {\omega }}+{\mathbf {\omega }}\times II.{\mathbf {\omega }}={\mathbf {M}}$

Tenembaum would accurately write the angular velocity vector as ${^{R}}{\mathbf {\omega }}^{A}$, however, to simplify the notation (and because we only have two reference frames here) we simply write ${\mathbf {\omega }}$ to refer to the angular velocity of the aircraft with respect to the earth. The matrix $II$ is the tensor of inertia, which for us is an identity matrix with $I_{xx}$, $I_{yy}$ and $I_{zz}$ in the main diagonal and 0's anywhere else.

In order to use the equations for numerical computations, it is necessary to write them in scalar form and to leave the derivatives alone at the left side:

Momentum

${\frac {^{A}d}{dt}}v_{x}={\frac {F_{x}}{m}}+g_{x}-\omega _{y}v_{z}+\omega _{z}v_{y}$

${\frac {^{A}d}{dt}}v_{y}={\frac {F_{y}}{m}}+g_{y}-\omega _{z}v_{x}+\omega _{x}v_{z}$

${\frac {^{A}d}{dt}}v_{z}={\frac {F_{z}}{m}}+g_{z}-\omega _{x}v_{y}+\omega _{y}v_{x}$

Angular momentum (Euler's equations)

${\frac {^{A}d}{dt}}\omega _{x}={\frac {1}{I_{xx}}}[M_{x}+(I_{yy}-I_{zz})\omega _{y}\omega _{z}]$

${\frac {^{A}d}{dt}}\omega _{y}={\frac {1}{I_{yy}}}[M_{y}+(I_{zz}-I_{xx})\omega _{z}\omega _{x}]$

${\frac {^{A}d}{dt}}\omega _{z}={\frac {1}{I_{zz}}}[M_{z}+(I_{xx}-I_{yy})\omega _{x}\omega _{y}]$

These six equations are not sufficient to describe the motion, since the gravity acceleration vector ${\mathbf {g}}$ also varies with time in the non inertial reference frame A. So we need to add the next three kinematic equations to the set:

${\frac {^{R}d}{dt}}{\mathbf {g}}={\mathbf {0}}\Rightarrow {\frac {^{A}d}{dt}}{\mathbf {g}}=-{\mathbf {\omega }}\times {\mathbf {g}}$

and thus, in scalar form:

Rotation of the gravity field

${\frac {^{A}d}{dt}}g_{x}=\omega _{z}g_{y}-\omega _{y}g_{z}$

${\frac {^{A}d}{dt}}g_{y}=\omega _{x}g_{z}-\omega _{z}g_{x}$

${\frac {^{A}d}{dt}}g_{z}=\omega _{y}g_{x}-\omega _{x}g_{y}$

Finally, if we want to output the position of the aircraft, then we need to integrate it in parallel also:

Position

${\frac {^{R}d}{dt}}P_{x}=v_{x}+\omega _{y}v_{z}-\omega _{z}v_{y}$

${\frac {^{R}d}{dt}}P_{y}=v_{y}+\omega _{z}v_{x}-\omega _{x}v_{z}$

${\frac {^{R}d}{dt}}P_{z}=v_{z}+\omega _{x}v_{y}-\omega _{y}v_{x}$

Note that this position will be relative to the aircraft reference coordinates, so to make it mean full, all points should be remapped to the global reference frame. The program does this when loading the results using the last orientation of the model.

Numerical integration

All the equations derived above form an initial value problem with 12 variables that can be solved using an appropriate numerical integration algorithm. One would probably think of the Runge-Kutta methods as first option, however they are not a good option since they require evaluating the forces at an intermediate time step, creating a complex scheme in our case (the UVLM method is quite heavy to be called continuously). Therefore, probably the best option is to use some kind of predictor corrector method that uses a fixed time step.

The numerical algorithm programmed in OpenVOGEL follows an scheme from Preidikman. This scheme is self-starting. It starts with the Euler method, then it switches to two Adams-Bashford and Adams-Moulton steps of increasing grade, and from the fourth step on it uses the Hamming method only.

In the next frames, vector X is the array (or better said, the record) containing all independent variables of the equations, and vector DX is the array containing the associated derivatives.

Time step 1

> Predict using initial derivatives (Euler)

 X(1) = X(0) + Dt * DX(0)

> Correct for K=0 to M (Modified Euler)

 Calculate new forces and derivatives...

 X(1) = X(0) + (0.5# * Dt) * (DX(0) + DX(1))

Time step 2

> Predict using previous derivatives (Adams-Bashford)

 X(2) = X(1) + (0.5# * Dt) * (3.0# * DX(1) - DX(0))

> Correct for K=0 to M (Adams-Moulton)

 Calculate new forces and derivatives...

 X(2) = X(1) + (Dt / 12.0#) * (5.0# * DX(2) + 8.0# * DX(1) - DX(0))

Time step 3

> Predict using previous derivatives (Adams-Bashford)

 X(3) = X(2) + (Dt / 12.0#) * (23.0# * DX(2) - 16.0# * DX(1) + 5.0# * DX(0))

> Correct for K=0 to M (Adams-Moulton)

 Calculate new forces and derivatives...

 X(3) = X(2) + Dt / 24.0# * (9.0# * DX(3) + 19.0# * DX(2) - 5.0# * DX(1) + DX(0))

 TE(3) = (9.0# / 121.0#) * (X(3) - XS)

Time steps S = 4 to N (Hamming)

> Predict using previous derivatives

 XP = X(S - 4) + (Dt * 4.0# / 3.0#) * (2.0 * DX(S - 1) - DX(S - 2) + 2.0 * DX(S - 3))
 
 X(S) = XP + 112.0# / 9.0# * TE(S - 1)

> Correct for K=0 to M

 Calculate new forces and derivatives...

 X(S) = (1.0# / 8.0#) * (9.0# * X(S - 1) - X(S - 3) + 3.0# * Dt * (DX(S) + 2.0# * DX(S - 1) - DX(S - 2)))
 
 TE(S) = (9.0# / 121.0#) * ((X(S) - XP))
 
 X(S) = X(S) - TE(S)

Validation

Validation of the numerical integration algorithm using a simple harmonic oscillator

Our VB.NET implementation of the algorithm has been validated for a simple 1-D harmonic oscillator. In such a system the force is a linear function of the position and the velocity:

$F_{x}=-K.P_{x}-C.v_{x}$

The analytic solution for the displacement considering only an initial velocity is :

$P_{x}(t)=e^{-\xi \omega t}v_{x}(0)sin(\omega _{d}t)$

where:

$\omega ={\sqrt {\frac {K}{M}}}$

$C_{c}=2{\sqrt {KM}}$

$\xi ={\frac {C}{C_{c}}}$

$\omega _{d}=\omega {\sqrt {1-\xi ^{2}}}$

The agreement of the numerical solution for this problem is excellent, as can be seen on the previous graphic. With 400 time steps and a time inteval of 0.025s, the simulation time was 10 seconds. Along the complete integration domain, the error could be maintained under 0.5%.

Validation cases

Validating software is an essential part of the development process. Software should not be use in practical applications without being confident that the results will not be far from reality. This chapter comprises the formal validation of OpenVOGEL's calculation core and is meant to provide sufficient evidence of the program's capability to predict aerodynamic loads, so as to support its application in real cases.

It is a fact that no software will ever provide an exact solution to a particular problem, and how intricate and complete our calculation model could be, there will always be some disagreement with reality. To mitigate this, the validation process will provide us an indication on the level of accuracy for different cases. For some models the accuracy might be very good, and for others it might be unacceptable. It could also be that the accuracy will be acceptable only under a certain subset of the problem's independent variables. Therefore it is important to study as many different configurations as possible so as to identify the situations that lay outside or at the boundary of the theoretical and numerical models.

The validation process can be split in two parts:

• Validation of the theoretical method: are we using the right assumptions?
• Validation of the numeric algorithms: is the code well written?

To answer these question, several test cases are going to be presented and analyzed in this section. We will start the validation with a wind tunnel test done in the year 1951 by NACA (currently NASA) at the Ames Aeronautical Laboratory. The model in this experiment is expected to be at the boundary of the our calculation capability since it treats a very low aspect ratio and high sweepback wing. For this case, the low subsonic speed might also give us an indication of how good a simple compressibility correction can be.

After having gained some confidence in the calculation core with the NACA model, we will go to the more conventional case of a low aspect ratio rectangular wing using data from the RAE-916 report (back to the year 1967).

As conclusion of both cases we will see how Tucan in combination with XFoil is very good at predicting the main aerodynamic forces on lifting surfaces up to moderate incidence angle.

The validation will then be focused on closed bodies. The first model will be a simple sphere, which will only be compared against theoretical results as Tucan does not handle flow separation. Finally, it is the idea to compare the static pressure prediction along a conventional fuselage using data from a more complex CFD model.

All validation cases are conform to calculation core version 2.1-2020.05 (this corresponds to the serial number prompted at the beginning of each calculation).

NACA RM-A51G31 technical report

Comparing results against the measurements in a wind tunnel is one of the best alternatives we have, because they were obtained using real air.NACA RM-A51G31 documents a very interesting wind tunnel test of a low aspect ratio wing at increasing Mach numbers. Since Tucan is intended for low speed aerodynamics, the lowest Mach of 0.25 has been selected for this test.

The model in Tucan consist of a flat lifting surface of the same overal dimensions as in the real model:

Top view of the Tucan model for NACA RM-A51G31 at 6°

Front view of the Tucan model for NACA RM-A51G31 at 6°

Property	Value	Value
Semi wingspan	46.67in	1.1854m
Root chord	41.47in	1.0533m
Tip chord	20.74in	0.5268m
Reference area	20.166ft²	1.8735m²
Dihedral	0°
Leading edge sweep-back	48.54°
Airfoil	NACA 64A010

Remarks

• In order to correct the compressibility effects, the Prandtl–Glauert transformation formula for 2D cases is used.

• The panels are uniformly spaced in span-wise direction and the number of panels is as follows:

Direction	Panels
Span-wise	44
Chord-wise	8
Total	352

• The simulation in Tucan was set up using the next parameters:

Property	Value
Wake length	60 panels
Steps	80
Step	0.015s
Airspeed	30m/s

Lift coefficient $C_{L}$

Lift coefficient versus incidence angle for the NACA test case. The wake is shed from the trailing edge (TE) and from the trailing edge plus wingtips (TE+WT).

In the next table the lift coefficient prediction from Tucan for the wake being shed from the trailing edge only (no wing tips) are compared against the NACA experiment. The values from the experiment have been extracted using pixels from the graph with the given accuracy.

	Tucan	Tucan	NACA
$\alpha$	M=0.00	M=0.25	M=0.25
0	0.0000	0.0000	0.0000
6	0.2969	0.3067	0.3047±0.005
12	0.5736	0.5924	0.6194±0.005

What we obtain as result is what we were expecting. The lift is predicted as a linear function of the incidence angle, but it does not drop after flow separation because the method we are using is unable to do this. The lift slope is accurately predicted and this gives the first indication that the algorithms are correct.

Drag coefficient $C_{D}$

Comparison of wind tunnel data against Tucan drag prediction (with and without wing tip wake shedding).

For the drag coefficient XFoil has been used to obtain the airfoil skin friction component at a Reynolds number of 4.000.000 to match the wind tunnel test conditions. The XFoil prediction at zero lift accounts for the total drag at zero incidence angle, and it is impressively accurate. As the incidence increases, the induced drag component predicted by Tucan takes over. The correlation is very good at low incidence, and fairly good at moderate incidence. At high incidence below the separation zone, the induced drag predicted by Tucan is considerably below the wind tunnel measurements and cannot be used.

It is interesting to note that for this model the drag prediction is much more consistent with reality when the wake shedding is done from the trailing edge and the wing tips (and not only from the trailing edge). This might have to do with the sharp sweep-back of the model (it is known that high sweep in delta wings causes leading edge flow separation), although it could also be that we are masking a different mechanisms by just adding more vorticity to the flow. This was briefly discussed in the tutorial about lifting surfaces: the extension of the wake shedding edge is an extra independent variable in the problem that we have to control manually (the so named "Kutta condition").

For the sake of practicty it does make sence to add more shedding if this seems to correlate better with results. In fact, scientists do this kind of things all the time when they don't have a solid evidence about something. Unfortunately, there is not enough data to know how the error is attributed to each of the two information sources (XFoil and Tucan's induced drag prediction algorithm). This could be done in the future with a CFD model. If the data extracted from the plot is correct, the laminar bubble predicted by XFoil seems to last too long, and maybe XFoil is thus also underpredicting the skin drag.

RAE-916 technical report

Technical report RAE-916 describes a series of wind tunnel experiments on very low aspect ratio wings. For this validation case the AF/1 wing case has been selected, which corresponds to a rectangular wing of aspect ratio 4 and symmetrical RAE-101 airfoil. The reported Reynolds number is 1.600.000, and for the Tucan model an XFoil polar at Reynolds 1.000.000 was used (expecting a slight disagreement there).

Property	Value	Value
Wingspan	96.00in	2.4384m
Root chord	24.00in	0.6096m
Tip chord	24.00in	0.6096m
Reference area	16.00ft²	1.4864m²
Dihedral	0°
Leading edge sweep-back	0.0°
Airfoil	RAE-101

Remarks

• For this model, no compressibility corrections were done (due to the low speed of 125ft/s).

• The panels are uniformly spaced in span-wise direction and the number of panels is as follows:

Direction	Panels
Span-wise	26
Chord-wise	8
Total	208

• The simulation in Tucan was set up using the next parameters:

Property	Value
Wake length	60 panels
Steps	80
Step	0.036s
Airspeed	38.0m/s

• The wing tips of the RAE model were a bit rounded, something that has not been modeled in Tucan (due to the meshing limitations).

Lift coefficient $C_{L}$

Lift coefficient prediction for the RAE-916 test case.

The predicted lift coefficient shows an excellent agreement with the experimental data. For this rectangular model the TE model also provides a more accurate description of the lift.

Drag coefficient $C_{D}$

Drag coefficient prediction for the RAE-916 test case.

For this rectangular wing the Tucan TE model provides a more accurate description of the drag than the TE+WT model (in contrast to the NACA model, where the wing tip wake shedding positively contributed to higher accuracy). This evidences that a full wake sheding is not always more realistic.

The two cases presented here suggest that the more the sweep-back the more wing tip shedding we need to improve the accuracy. However this cannot be taken as conclusion without analyzing more cases in between and higher aspect ratios.

Thick bodies: potential flow around a sphere

The potential flow around a sphere can be solved analytically, and it is therefore a good reference for testing thick bodies. The analytical solution is

$C_{p}=1-9/4cos(\theta )$.

Therefore, we should expect 1 at the stagnation point, and -1.25 at the top.

The current algorithm used to find the local surface velocity is based in a estimation of the circulation gradient using least squares. Based on that velocity, the pressure is computed using Bernoulli's equation. The next graph shows how Tucan is currently able of approaching the Cp over the surface. The model consisted of 2400 quadrilateral panels. These results correspond to kernel versions 2.0 and higher (before that, the algorithm was not suitable).