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

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

尾迹寿命和特殊建模特性

计算核心的一个显着特征是，每个尾迹都可以采用自己的寿命。这允许用户控制每个尾迹的长度以及边界条件。例如，对于正常尾迹，在机翼后缘与机身合并的确切位置，机翼根部会出现问题。从那里脱落尾迹节点不是一个好主意，因为在表面附近引入了虚假速度分量。尽管如此，我们仍然可以通过引入一个 0 CuttingStep 的尾迹，在这个后缘的小部分区域强制执行库塔条件，而无需任何额外代码。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 标量来计算。中间点始终使用面板节点的平均坐标来计算。

涡流

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

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

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

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

全局邻接调查

为了计算细长涡环上的压力跃变，需要通过向量求和侧涡来计算局部涡量向量。此操作需要知道哪些涡环与每个涡环的每一侧相邻。负责此全局邻接调查的过程称为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}$ 指的是通过使用从原始脱落边缘继承的涡旋强度，利用毕奥-萨伐尔定律获得的第 $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}}]$