# The calculation core

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

## 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 ${\displaystyle \mathbf {u} ,\mathbf {v} ,\mathbf {w} }$ is created. For quadrilateral panels, ${\displaystyle \mathbf {u} }$ points in the direction of the first diagonal (from node 1 to node 3), ${\displaystyle \mathbf {w} }$ is normal to both diagonals, and ${\displaystyle \mathbf {v} }$ is normal to the other two vectors. For triangular panels, ${\displaystyle \mathbf {u} }$ points in the direction of the first edge (from node 1 to node 2), ${\displaystyle \mathbf {w} }$ is normal to the surface (which is always flat), and ${\displaystyle \mathbf {v} }$ is normal to the other two vectors. In all cases, vector ${\displaystyle \mathbf {v} }$ 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

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.

#### Vortices and the slender surface pressure jump

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.

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

${\displaystyle \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:

${\displaystyle \mathbf {V} (\mathbf {R} ,t)={\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. ${\displaystyle k=3}$ for triangles, and ${\displaystyle k=4}$ for quadrilaterals). The sub-indices (i,j) refer to the first and second nodes of the side ${\displaystyle 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 ${\displaystyle \mathbf {L} _{ji}}$ are the side segments given by ${\displaystyle \mathbf {L} _{ji}=\mathbf {R} _{j}-\mathbf {R} _{i}}$, where ${\displaystyle \mathbf {R} _{i}}$ and ${\displaystyle \mathbf {R} _{j}}$ are the position vectors of the two nodes on side ${\displaystyle k}$. Vector ${\displaystyle \mathbf {r} _{i}}$ is the relative position of the targeting point ${\displaystyle \mathbf {R} }$ respect to ${\displaystyle \mathbf {R} _{i}}$, that is to say, ${\displaystyle \mathbf {r} _{i}=\mathbf {R} -\mathbf {R} _{i}}$, and the unit vectors ${\displaystyle \mathbf {e} _{i}}$ and ${\displaystyle \mathbf {e} _{j}}$ point from each respective side node to the targeting point.

### Doublet panels: potential

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

${\displaystyle \phi (\mathbf {R} ,t)={\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. ${\displaystyle k=3}$ for triangles, and ${\displaystyle k=4}$ for quadrilaterals). The sub-indices (i,j) refer to the first and second nodes of the side ${\displaystyle 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:

${\displaystyle m_{ij}={\frac {y_{j}-y_{i}}{x_{j}-x_{i}}}}$

${\displaystyle r_{i}={\sqrt {(x-x_{i})^{2}+(y-y_{i})^{2}+z^{2}}}}$

${\displaystyle e_{i}=(x-x_{i})^{2}+z^{2}}$

### 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 [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:

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

### Sink/source panels: velocity

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

### Dirichlet boundary conditions

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:

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

In this expression, ${\displaystyle \phi _{G}}$ is the potential associated to the bounded source panels, ${\displaystyle \phi _{S}}$ is the potential associated to the bounded sink/source panels and ${\displaystyle \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 ${\displaystyle \mathbf {M} _{S}}$. This matrix can then be reused to compute the total sink/source potential when necessary as follows:

${\displaystyle \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 ${\displaystyle \mathbf {G} }$. Similarly to the sources, we can write a matrix of unit-intensity potential ${\displaystyle \mathbf {M} _{G}}$ so that:

${\displaystyle \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:

${\displaystyle \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

${\displaystyle \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.

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

${\displaystyle 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}$

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:

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

The doublets cross-flow becomes:

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

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

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

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

1. "Low speed aerodynamics", Allen Plotkin, Joseph Katz
2. "Low speed aerodynamics", Allen Plotkin, Joseph Katz
3. "Low speed aerodynamics", Allen Plotkin, Joseph Katz