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