The aerodynamic model in OpenVOGEL

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

The aerodynamic calculation core in classes

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Wake lifetime and special modeling features

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

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

Building blocks: vortices, doublets and sources

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

Public Interface VortexRing
   [...]

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

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

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

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

   [...]
End Interface

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

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

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

Local coordinates

Projection of a quad panel into a plane.

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

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

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

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

Collocation points

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

Vortices

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

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

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

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

Global adjacency survey

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

Simulation of an extended flap by a second surface.

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

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

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

The math problem

Introduction

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

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

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

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

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

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

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

Doublet panels: velocity

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

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

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

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

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

Doublet panels: potential

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

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

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

The formulation is completed with the next definitions:

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

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

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

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

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

Sink/source panels: potential

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

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

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

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

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

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

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

Sink/source panels: velocity

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

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

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

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

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

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

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

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

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

Dirichlet boundary conditions

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Local surface velocity in thick closed bodies

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

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

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

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

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

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

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

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

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

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

4. Finally we will have:

$dG/du=A_{1}$

$dG/dv=A_{2}$

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

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

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

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

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

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

Pressure coefficient in thick bodies

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

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

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

Neumann boundary conditions

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

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

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

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

with each component as follows:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Mixed boundary conditions

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

Reference

↑ "Low speed aerodynamics", Allen Plotkin, Joseph Katz
↑ "Low speed aerodynamics", Allen Plotkin, Joseph Katz
↑ "Low speed aerodynamics", Allen Plotkin, Joseph Katz

[1] "Low speed aerodynamics", Allen Plotkin, Joseph Katz

[2] "Low speed aerodynamics", Allen Plotkin, Joseph Katz

[3] "Low speed aerodynamics", Allen Plotkin, Joseph Katz

[1]

[2]

[3]