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