OpenVOGEL中的空气动力学模型

计算核心（CC）是程序中处理计算算法的部分。正如关于源代码的章节中所解释的，OpenVOGEL的编写方式使得计算核心独立于可视化模型。计算模型与可视化模型唯一的相遇时刻是在计算启动时，此时后者被转换为前者。计算核心也纯粹基于面向对象，而面板法非常适合此方法。事实上，每个面板类型都是通过一个类来实现的。类似地，用于气动弹性分析结构部分的有限元方法也以面向对象的方式进行编程。我们将从本章开始，看看空气动力学模型的主要组成部分。之后，我们将更详细地研究空气动力学问题的数学原理。

类中的空气动力学计算核心

OpenVOGEL通过引入两种主要基本类型来构建面板法：Vortex 类和VotexRing 接口。Vortex 类代表一条直线涡流段，而VortexRing 接口代表扁平面板。后者由两个类实现：VortexRing3，代表三角形，和VortexRing4，代表四边形。所有这些元素类型都基于重要的节点概念，该概念由Node 类表示。节点实际上是在流体域中跟踪的唯一物质点。

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

在对象的分类中，我们发现格架，它由通用Lattice 类表示。格架是一个由VortexRing 和/或Vortex 元素组成的集合，这些元素连接了Node 对象的云。

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

这意味着一个格架可以包含三角形和四边形面板的混合物，加上一个相互连接的涡流网络。格架公开方法来计算其组件在空间中给定点处的诱导速度，这些方法用于构建空气动力学影响系数和脱落尾流。

代码中还有两个从Lattice 类派生的额外类。它们是Wake 和BoundedLattice 类，它们具有非常特定的含义。Wake 是一个代表自由涡流片的格架，而BoundedLattice 是一个代表固体表面（或更确切地说，是薄边界层）的格架。

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

边界格架包含一个尾流列表，而尾流包含对从边界格架上脱落到流体域中的原始边缘面板的引用。这些原始面板提供了循环的值，该值将表征尾流环或涡流在其整个生命周期中的特征。要声明一个原始值，您必须引入脱落边缘处的节点索引，然后按顺序引入相应面板的索引。

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

下图清楚地显示了原始边缘的定义。从这张图中也可以清楚地看出，原始边缘并不局限于机翼的根部和翼尖，而是可以在任何地方定义，包括前缘。

尾流寿命和特殊建模特征

计算核心的一个显著特征是，每个尾流都可以采用自己的寿命。这允许用户控制每个尾流的长度以及边界条件。例如，对于普通尾流，在机翼后缘与机身合并的机翼根部会出现一个问题。从那里脱落尾流节点不是一个好主意，因为会在表面附近引入虚假的速度分量。尽管如此，我们仍然可以通过引入一个CuttingStep 为 0 的尾流，在这个后缘的小部分区域强制实施库塔条件，而无需任何额外的代码。OpenVOGEL在转换模块中自动执行此操作，而内核并不知道该尾流的实际用途。

另一个需要内核意识到的问题是内循环的抑制。这个问题再次出现在机翼根部，但方向是来流方向。靠近机身的涡流在现实中并不能自由移动，因为它们非常靠近机身。因此，为了避免尾流的内侧部分滚动，流向机翼根部的第一条尾流涡线必须被抑制。内核可以通过激活SupressInnerCircuation 属性来为每个尾流执行此操作。

构建块：涡流、偶极子和源

现在我们已经介绍了空气动力学计算核心的主要类，我们可以更深入地了解。但在开始讨论主要计算算法之前，我们可能需要了解涡流和涡流环的功能。如前所述，VortexRing 派生类以两种形式存在，分别具有三个节点和四个节点。但这两种类中的任何一个也可以以两种不同的方式运行：作为扁平的恒定偶极子或作为恒定源分布区域。为此，这两个类实现了一些方法来计算偶极子速度、偶极子势、源速度和源势，这些值在空间中给定点处。VortexRing 派生类中封装的下一组函数构成了计算核心的基本资源，因为它们允许确定空气动力学影响系数。

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

当一个面板表现为恒定偶极子分布时，则使用前两种方法。速度影响对于计算空间中任何给定点处的诱导速度很有用。如果参数WithG 设置为 true，则使用局部循环值，结果是绝对速度。如果WithG 设置为 false，则使用单位循环，结果是单位循环的诱导速度。在相同的行为下，第二个函数返回空间中任何给定点处的势函数值，而参数WithG 的作用与之前相同。当面板表现为恒定源分布时，第三个函数返回空间中任何点处的势函数值。参数WithS 与WithG 类似，指示是否必须使用单位源/汇强度。

你现在可能想知道，在什么情况下我们使用面板的双极子行为或源行为，以及什么时候我们需要计算诱导速度或局部势。答案是，这取决于所施加的边界条件。在细长表面面板上，我们总是施加 Neumann 边界条件，而在厚表面面板上，我们总是施加 Dirichlet 边界条件。因此，整个边界条件问题自然地被分为两部分。

当 Neumann 边界条件施加在面板上时，必须在位于表面的中间控制点计算诱导速度。因此，我们将扫描所有表面和尾流面板，并利用它们的双极子和源行为来计算局部诱导速度。另一方面，当 Dirichlet 边界条件施加在面板上时，必须在位于表面正下方的内部控制点计算势，为此，我们扫描所有表面和尾流面板来计算速度势。

局部坐标

在 OpenVOGEL 中，VortexRing 派生类在局部坐标系中实现势函数的计算。为此，创建了一个正交向量 $\mathbf {u} ,\mathbf {v} ,\mathbf {w}$ 的参考基底。对于四边形面板， $\mathbf {u}$ 指向第一个对角线（从节点 1 到节点 3）的方向， $\mathbf {w}$ 与两个对角线都垂直，而 $\mathbf {v}$ 与另外两个向量垂直。对于三角形面板， $\mathbf {u}$ 指向第一条边（从节点 1 到节点 2）的方向， $\mathbf {w}$ 与表面垂直（总是平坦的），而 $\mathbf {v}$ 与另外两个向量垂直。在所有情况下，向量 $\mathbf {v}$ 总是选择为使得基底为右手系。当一个面板被声明为反转时，所有向量都只是翻转。

在接下来的章节中，我们将使用符号 $(u,v,w)$ 来表示在这个局部基底中相对于面板控制点的坐标。大多数参考文献使用 $(x,y,z)$ 来表示，但使用与全局坐标系相同的符号会造成混淆（指的是方向向量 $\mathbf {i} ,\mathbf {j} ,\mathbf {k}$ ）。当我们指的是速度分量时，我们将使用 $(V_{x},V_{y},V_{z})$ 或 $(V_{u},V_{v},V_{w})$ 。然而，在代码中，我们被迫使用字母 x、y 和 z，因为我们在所有情况下都使用相同的 Vector3 类。

使用这个局部基底，顶点的位置以二维方式重新计算并缓存以提高性能，以便在执行过程中不必连续重新计算它们（例如，在 Katz & Plotkin 的 FORTRAN 示例中就是这样）。出于同样的原因，边的长度也缓存。

请注意，局部基底的声明不是 VortexRing 接口的必要条件。实际上，接口的任何实现都可以在其私有部分自由选择自己的投影方法。在当前版本的 CC 中，局部基底方法也已在三角形中实现，无论它们是否始终是平坦的。

配置点

如前所述，边界条件问题需要声明多个配置点。在 OpenVOGEL 中，每个面板计算三个配置点：内部控制点、中间控制点和外部控制点。内部点仅用于计算厚体的内部势。中间点用作局部坐标的原点（以完成投影面板的定义），并用于计算细长表面的速度。最后，外部点用于评估厚体上的速度。内部点和外部点的位置使用中间点和法线方向乘以 epsilon 标量来计算。中间点始终使用面板节点的平均坐标计算。

涡流

虽然整个计算核心基本上是基于VortexRings（面板）的，但有时使用更基本类型（Vortex）是有原因的。这种情况出现在脱落尾流和需要施加 Neumann 边界条件时，也就是说，当模型中存在细长表面时。在这种情况下，需要在细长面板的控制点和尾流的节点处计算感应速度，为此使用Vortices 更有效率。事实上，当我们使用VortexRing 在某处计算速度时，我们会绕着它的三个或四个边界段循环。如果我们对两个相邻的环这样做，那么我们会访问它们共有的边缘段两次。但是，通过实现涡流，我们可以避免这种情况，并将计算次数减少到几乎一半。因此，我们使用Vortex 元素的主要原因是计算效率。

利用这一点的最佳方法是为网格配备两个并行结构，一个VortexRings 网格和一个Vortices 网格，并在最方便的情况下使用其中一个。在 OpenVOGEL 中，计算核心将始终使用VortexRings 和边界网格中的Vortices 的并行网格。在Wakes 上，情况有所不同：它将始终使用Vortices 网格，并且VortextRings 仅在需要在某处应用 Dirichlet 边界条件时才会使用（因为它们是计算势所需的）。通过这样做，我们获得了双重收益：我们通过使用VortexRings 减少了矩阵问题的大小，因为VortexRings 的数量总是少于Vortices，并且我们通过使用Vortices 减少了构建右侧和脱落尾流所需的计算次数，因为它们需要更少的运算次数。

当我们使用涡流的并行网格时，为每个Vortex 提供对两个或三个相邻VortexRings 的引用，以及有关其方向的信息（说明如何解释其环流的符号），非常有效。有了所有这些数据，每个涡流的强度可以通过简单求和来评估。基于VortexRings 生成网格的涡流，并搜索和分配每个Vortex 的相邻VortexRings 的过程称为PopulateVortices。此方法位于Lattice 类中（因此对所有网格都是通用的）。请注意，该过程不需要扫描所有网格，但必须能够找到至少 3 个相邻环（因为这种情况发生在机翼-机身锚点处）。

对于尾流，情况再次不同。由于尾流环和涡流的环流在其整个生命周期中保持不变，因此它们的环流在尾流脱落过程中的创建时直接分配。

全局邻接调查

为了计算细长VortexRing 上的压力跳跃，需要通过向量求和侧涡来计算局部涡量向量。此操作需要知道哪些VortexRings 与每个VortexRing 的每一侧相邻。负责此全局邻接调查的过程称为FindSurroundingRingsGlobally，它位于源文件AeroTools.CalculationModel.Solver_Calculations 中。此过程必须扫描所有网格，因为在某些接口（例如机翼-机身锚点）处，环可能不共享相同的节点。当两个节点之间的距离小于给定容差时，它们被认为位于同一点。

全局邻接调查过程保证即使面板不共享相同的节点，相邻面板上的压力跳跃也能被正确计算。唯一的要求是连接共享边的节点彼此足够接近（比Settings 中声明的SurveyTolerance 参数更近）。此功能的实际应用是弗拉普襟翼 或控制表面的建模。如果此过程不包含在计算核心中，则两个表面沿其环流的连续性将被破坏，并且每个面板都会在共享边上看到一个开口。下图代表了一个实际示例。请注意脱落边缘是如何声明的，以及升力是如何在襟翼铰链线后面增加的。

除了所有这些之外，全局邻接调查过程对于使用环流梯度计算厚体上的局部速度也是必需的（这将在后面关于 Dirichlet 边界条件的部分中进行解释）。

数学问题

介绍

由于 OpenVOGEL 处理多种类型的面板和边界条件，因此数学和实现的算法比仅处理涡流环所需的算法复杂得多。

处理面板时的基本问题是找出导致其表面没有法向流动的环流。对于细长面板，此要求通过说明流体在控制点处的法向分量必须等于零来满足。这就是所谓的Neumann 边界条件。由于与特定点处的涡流环相关的速度线性取决于该环的环流，因此可以写出线性方程组，以便在求解时，它将提供取消每个控制点处法向速度的环流值。

现在，当需要对诸如机身之类的封闭物体进行建模时，以前的方法有时会产生一个病态系统，无法求解。解决方法是包含汇/源面板，并通过说明物体内部的总势必须等于给定的常数值（通常选择为零）来施加非穿透条件。这就是所谓的Dirichlet 边界条件。请注意，这是一种非常不同的问题类型，需要计算与偶极子分布和汇/源分布相关的势。此外，为了评估周围的流速场，我们还必须事后计算这些汇/源面板的速度影响（局部表面速度的处理方式不同，如下所述）。因此，在最一般的情况下，需要四个基本算法

计算涡流环的速度
计算偶极子平面分布的势
计算汇/源平面分布的速度
计算汇/源平面分布的势

因此，在我们能够制定完整的边界条件并求解面板的环流之前，有必要写出与面板相关的空间中任意点的势和势导数（速度）的数学表达式，无论是均匀汇/源分布还是均匀偶极子分布。

请注意，我们在势流下，因此下一个关系成立

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

偶极子面板：速度

偶极子面板代表一个由三个或四个直涡段组成的涡环。面板对空间中任何点速度的贡献由著名的Biot-Savart 公式给出。积分方程得出

$\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})]}$

求和扩展到面板的每一侧（例如， $k=3$ 对于三角形，以及 $k=4$ 对于四边形）。子索引 (i,j) 指的是边 $k$ 的第一个和第二个节点，因此它们采用值 $(1,2);(2,3);(3,1)$ 对于三角形和 $(1,2);(2,3);(3,4);(4,1)$ 对于四边形。向量 $\mathbf {L} _{ji}$ 是由 $\mathbf {L} _{ji}=\mathbf {R} _{j}-\mathbf {R} _{i}$ 给出的边段，其中 $\mathbf {R} _{i}$ 和 $\mathbf {R} _{j}$ 是边 $k$ 上两个节点的位置向量。向量 $\mathbf {r} _{i}$ 是目标点 $\mathbf {R}$ 相对于 $\mathbf {R} _{i}$ 的相对位置，也就是说， $\mathbf {r} _{i}=\mathbf {R} -\mathbf {R} _{i}$ ，而单位向量 $\mathbf {e} _{i}$ 和 $\mathbf {e} _{j}$ 分别从每个边节点指向目标点。

涡旋段速度的 VB.NET 数值实现如下。

''' <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

双层面板：势

不深入细节，由于均匀的双层分布，平板的势可以从 ^[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}}})]$

该公式用局部坐标表示，求和扩展到面板的每一侧（例如，对于三角形， $k=3$ ，对于四边形， $k=4$ ）。下标(i,j)指代边 $k$ 的第一个和第二个节点，因此它们取值 $(1,2);(2,3);(3,1)$ （三角形）和 $(1,2);(2,3);(3,4);(4,1)$ （四边形）。请注意，上述等式在数值算法中被进一步简化，以避免计算两个反正切函数并提高性能。

该公式通过以下定义完成

$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}$

上述等式的 VB.NET 数值实现如下。

''' <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

汇/源面板：势

计算与恒定汇/源面板和涡环（恒定双极子）在空间中给定点相关的速度势非常复杂。OpenVOGEL 通过将四边形面板投影到其法线方向（根据局部基底）来降低复杂度级别，因此它们可以由一个平面面板表示（三角形面板不需要投影，因为它们始终是平面的，尽管它们的处理方式类似）。在这种情况下，速度势的计算可以在该参考文献中找到^[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}}})]$

当然，这种方法也存在一些泄漏问题，这取决于面板及其平面对应物之间的匹配程度。因此，重要的是要提供一个网格，使四边形面板不会过分扭曲。

在使用汇/源面板时，与其相关的强度不是未知量。事实上，当目标内部势能为零时，汇/源强度只需设置为自由气流速度的法线分量。

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

三角形面板势能的 VB.NET 数值实现如下。

''' <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

汇/源面板：速度

与汇/源面板相关的速度可以在参考文献 ^[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}}})]$

全局坐标系中的速度可以通过将三个正交投影进行求和来构建。

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

源速度影响的 VB.NET 数值实现如下。

''' <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

请注意，在代码的最后，局部速度分量（在局部坐标系中）使用局部面板正交基投影回全局坐标系。

狄利克雷边界条件

所谓 *狄利克雷边界条件* 施加在厚闭合物体的表面，要求物体表面下的势能保持恒定。为了用数学术语写出这一点，我们将与有界偶极子和源相关的势能分解如下

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

在这个表达式中， $\phi _{G}$ 是与有界偶极子面板相关的势能， $\phi _{S}$ 是与有界汇/源面板相关的势能， $\phi _{W}$ 是与尾迹相关的势能（根据定义，尾迹只包含偶极子）。将常数设置为零的偏好与之前的假设有关，即汇/源强度由自由流速度的法向分量决定。

由于每个有界面板的汇/源强度是预先知道的，对于刚性面板（即连接到不移动节点的面板），可以在每个有界内部控制点上计算一次单位势能并存储在汇/源势能矩阵 $\mathbf {M} _{S}$ 中。然后，可以根据需要重复使用此矩阵来计算总的汇/源势能，如下所示

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

请注意，这种技术可能不太适合变形格子。但是，在 OpenVOGEL 中，不会遇到这个问题，因为气动弹性问题仅保留给细长面板。因此，无论何时我们使用厚物体，相应面板的形状在整个计算过程中始终保持不变。

狄利克雷问题的未知数是包围物体的每个面板的偶极子强度，这里用向量 $\mathbf {G}$ 表示。类似于源，我们可以写一个单位强度势能的矩阵 $\mathbf {M} _{G}$ ，以便

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

最后，需要加上与尾迹相关的势能。由于尾迹不断改变形状，我们更倾向于将它们的势能作为单个贡献的总和，而不是构建矩阵

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

由于尾迹是从升力面上脱落的，它们的势能只与已知强度的恒定偶极子面板有关。当尾迹面板距离物体足够远时，通常的做法是通过点偶极子来近似它们对势能的贡献，点偶极子具有更简单的数学表达式。此解称为 *远场势能*。但是，当前版本的 CC 尚未包含此功能。始终使用完整的方程。

然后，狄利克雷问题由以下线性方程组表示

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

这需要与诺伊曼问题（见下一节）一起在全局矩阵问题中进行组装和求解。

封闭厚物体表面局部速度

一旦找到环流 $\mathbf {G}$ ，局部速度不是通过将速度影响相加来计算的。实际应用表明，这种方法不能得到正确的结果，因为它无法解决由于表面环流引起的局部切向速度。相反，局部速度是通过局部二维环流梯度来求解的。一般来说，我们有

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

这里的问题是如何确定局部导数 $dG/du$ 和 $dG/dv$ 。当然，对于函数 $G(u,v)$ ，不存在解析公式，因为假定每个面板上的环流是恒定的。因此，我们需要使用相邻面板的信息。

Tucan 通过使用最小二乘法，利用相邻面板的环流拟合微分函数，以相对简单的方式解决了这个问题。该方法的灵感来自中心差分法（以及中值定理），但这种方法的优点是，只要存在至少两个或多个非共线的相邻面板（对于封闭物体，这种情况应该总是存在的），它就可以应用于任何类型的面板，而无需任何修改。该算法为：

1. 利用三个或四个相邻面板以及面板本身，构建矩阵 $\mathbf {M}$ ，如下所示

$\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]$

其中 $n$ 是样本点的总数（通常对于三角形和四边形面板分别为 4 或 5 个）。

2. 构建右手边向量 $\mathbf {B}$ ，如下所示

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

3. 求解 3x3 线性方程组 $\mathbf {M} \mathbf {A} =\mathbf {B}$ .

4. 最后，我们将有

$dG/du=A_{1}$

$dG/dv=A_{2}$

解向量中的第三个分量是平均循环值，可以丢弃。

接下来的 VB.NET 代码展示了该算法的实现方式

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

需要注意的是，当一个面板在一侧有多于两个相邻面板时，该侧将被忽略，因为假设其中一个面板正在阻挡另一个面板（这与作为锚点连接到主体的细长面板的情况一致）。该方法的另一个要求是，必须找到模型中的所有相邻面板（使用之前描述的全局邻接调查过程）。

上述数值方法实际上产生了所需的速率，这听起来可能很奇怪。但是，正如之前解释的那样，技巧在于使用相邻面板的循环来最佳拟合以目标面板为中心的二维微分（即线性）函数。该算法的局限性可能在于表面的曲率。需要注意的是，使用每个相对位置向量在每个局部方向上的投影来逼近梯度，因此它采取了到相邻控制点的捷径，而不是使用实际表面坐标。

然而，在球体上验证算法（球体已知解析解）表明，得到的局部速度和 $C_{P}$ 的预测精度非常高（请参阅下图）。

厚物体中的压力系数

确定局部表面速度后，可以使用伯努利方程计算压力系数。对于稳定情况，我们有

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

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

诺伊曼边界条件

在细长表面上，我们应用所谓的诺伊曼边界条件，这些条件通过说明穿过表面的法向速度分量必须为零（即没有气流穿过表面）来表达。然后，我们可以将控制点处的总速度分解为不同的分量，并将其投影到面板法向量上以获得净横流

$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$

此方程必须针对每个细长面板重复。我们将 $N_{G}$ 命名为细长面板的总数，这意味着将存在一个 $N_{G}\times N_{G}$ 方程组。与狄利克雷边界条件的处理方式类似，这里我们可以为双极子构建一个影响横流矩阵，为源构建另一个矩阵，并将尾迹横流在求和中分离。如前几章所述，与尾迹相关的速度 $\mathbf {V} _{W}(\mathbf {R} _{i})$ 可以通过使用 $N_{WV}$ 尾迹涡流更好地计算。我们可以将尾迹相关的横流向量写成以下形式

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

每个分量如下

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

这里 $\mathbf {V} _{Vk}$ 代表由 Biot-Savart 方程得到的第 $k$ 个尾迹涡的诱导速度，该速度使用从原始脱落边继承的涡强度计算得出。

自由来流的横流分量可以写成以下形式：

$\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}}]$