# Open VOGEL/Kernel

## Contents

# The calculation core[edit]

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[edit]

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
[...]
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.

## Building blocks: vortices, doublets and sources[edit]

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[edit]

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 is created. For quadrilateral panels, points in the direction of the first diagonal (from node 1 to node 3), is normal to both diagonals, and is normal to the other two vectors. For triangular panels, points in the direction of the first edge (from node 1 to node 2), is normal to the surface (which is always flat), and is normal to the other two vectors. In all cases, vector is always chosen so that the basis is dextrorotation. When a panels is declared as *reversed*, all vectors are just flipped.

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[edit]

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[edit]

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.

#### Vortices and the slender surface pressure jump[edit]

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.

## The math problem[edit]

### Introduction[edit]

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 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 local pressure coefficient the velocity induced by the sink/source distributions will also have to be calculated afterwards. 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:

### Doublet panels: velocity[edit]

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:

The summation extends to each side of the panel (i.e. for triangles, and for quadrilaterals). The sub-indices *(i,j)* refer to the first and second nodes of the side , 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 are the side segments given by , where and are the position vectors of the two nodes on side . Vector is the relative position of the targeting point respect to , that is to say, , and the unit vectors and point from each respective side node to the targeting point.

### Doublet panels: potential[edit]

Without going into details, the potential of a flat panel due to a uniform doublet distribution can by obtained from ^{[1]}:

This formulation is written in local coordinates and the summation extends to each side of the panel (i.e. for triangles, and for quadrilaterals). The sub-indices *(i,j)* refer to the first and second nodes of the side , 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:

### Sink/source panels: potential[edit]

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 ^{[2]}. 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:

### Sink/source panels: velocity[edit]

The velocity associated to a sink/source panel can be found in reference ^{[3]}.

### Dirichlet boundary conditions[edit]

The Dirichlet boundary conditions are imposed on the surface of thick closed bodies, and they can be expressed by stating that the potential immediately under the surface must be constant. If we separate the potential associated to bounded doublets and sources we can write:

In this expression, is the potential associated to the bounded source panels, is the potential associated to the bounded sink/source panels and 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 . This matrix can then be reused to compute the total sink/source potential when necessary as follows:

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 . Similarly to the sources, we can write a matrix of unit-intensity potential so that:

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:

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

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

### Neumann boundary conditions[edit]

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):

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.

The sources cross-flow becomes:

The doublets cross-flow becomes:

The unknown of the problem are the source intensities associated to the slender panels (vector ), which can be found after solving the next system of linear equations:

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