# OpenVOGEL/Aerodynamics

# The aerodynamic model in OpenVOGEL[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
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.

### Wake lifetime and special modeling features[edit]

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

In the following sections we will be using the notation to refer to coordinates in this local basis and respect to the panel control point. Most of the bibliography uses for this, but it is confusing to use the same notation as for the global coordinates system (refering to direction vectors ). When we refer to velocity components, we will use or . 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[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.

### Global adjacency survey[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.

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

### 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 for triangles and 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.

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[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 for triangles and 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:

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

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

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

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

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 referece 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[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 doublet 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.

#### Local surface velocity in thick closed bodies[edit]

Once the circulations 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:

The problem here is how to determine the local derivatives and . Of course there is no analytic formula for the function , 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 as follows:

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

2. Build right hand side vector as follows:

3. Solve the 3x3 system of linear equations .

4. Finally we will have:

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 are predicted with very good accuracy (see picture here next).

#### Pressure coefficient in thick bodies[edit]

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

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

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

This equation must be repeated for each slender panel. We will name to the total number of slender panels, meaning that there will be a system of 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 can better be calculated using the wake vortices. We can write a wake related cross flow vector as follows:

with each component as follows:

By we mean here the influence velocity of the wake vortex 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:

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

The matrix 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 and each component can be written as follows:

By we mean here the influence velocity associated to the flat sink-source distribution over panel , which can be obtained from the formulas that were written in the previous paragraphs.

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

The matrix contains the influence cross flow of the slender surface panels due to the panels dipole strength (circulation). Therefore, the size of this matrix is and each component can be written as follows:

By we mean here the influence velocity associated to the flat dipole distribution over panel (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 ), 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.

### Mixed boundary conditions[edit]

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.