Make LiftingLine and StripwiseELement inherit from AbstractBody and AbstractElement; postprocess linear solution

Author: EdoAlvarezR

src/FLOWPanel.jl

@@ -64,7 +64,7 @@ for header_name in ["elements", "linearsolver",
                     "abstractbody", "nonliftingbody",
                     "abstractliftingbody", "liftingbody",
                     "multibody",
-                    "liftingline_stripwise", "liftingline",
+                    "liftingline_stripwise", "liftingline", "liftingline_postprocess",
                     "utils", "postprocess",
                     # "fmm"
                     ]

src/FLOWPanel_liftingline.jl

@@ -15,16 +15,17 @@ import LaTeXStrings: @L_str
 # LIFTING LINE STRUCT
 ################################################################################
 struct LiftingLine{ R<:Number, 
-                    S<:StripwiseElement, 
+                    S<:StripwiseElement, N,
                     VectorType<:AbstractVector{R}, 
                     MatrixType<:AbstractMatrix{R}, 
                     TensorType<:AbstractArray{R, 3},
-                    LI<:LinearIndices}
+                    LI<:LinearIndices} <: AbstractBody{S, N}
 
     # Internal properties
     grid::gt.Grid                               # Flat-geometry grid
     linearindices::LI                           # Linear indices of grid.nodes where linearindices[i, j] 
                                                 # is TE (j==1) or LE (j==2) of the i-th row of nodes
+    fields::Vector{String}                      # Available fields (solutions)
 
     ypositions::Vector{Float64}                 # Non-dimensional y-position of nodes, 2*y/b
 
@@ -44,12 +45,13 @@ struct LiftingLine{ R<:Number,
                                                 # j-th semi-infinite filament (1==a, 2==b) of the
                                                 # n-th horeseshoe
 
-    midpoints::MatrixType                       # Midpoint along bound vortex for probing the velocity
+    midpoints::MatrixType                       # Midpoint along lifting line for probing the velocity
     controlpoints::MatrixType                   # Control point of each horseshoe
     normals::MatrixType                         # Normal of each horseshoe
 
     aoas::VectorType                            # (deg) angle of attack seen by each stripwise element
     Gammas::VectorType                          # Circulation of each horseshoe
+    Us::MatrixType                              # Velocity at each midpoint
 
     G::MatrixType                               # Geometry matrix in linear system of equations
     RHS::VectorType                             # Right-hand-side vector in linear system of equations
@@ -116,11 +118,13 @@ struct LiftingLine{ R<:Number,
         normals = MatrixType(undef, 3, nelements)
         aoas = VectorType(undef, nelements)
         Gammas = VectorType(undef, nelements)
+        Us = MatrixType(undef, 3, nelements)
         G = MatrixType(undef, nelements, nelements)
         RHS = VectorType(undef, nelements)
 
         aoas .= 0
         Gammas .= 0
+        Us .= 0
         G .= 0
         RHS .= 0
 
@@ -139,17 +143,17 @@ struct LiftingLine{ R<:Number,
 
         calc_Dinfs!(Dinfs, initial_Uinf, nelements)
 
-
+        S = eltype(elements)
         new{R,
-            eltype(elements),
+            S, _count(S),
             VectorType, MatrixType, TensorType, 
             typeof(linearindices)}(
-                                grid, linearindices,
+                                grid, linearindices, String[],
                                 ypositions, 
                                 nelements, elements,
                                 aerocenters, horseshoes, Dinfs, 
                                 midpoints, controlpoints, normals,
-                                aoas, Gammas,
+                                aoas, Gammas, Us,
                                 G, RHS,
                                 kerneloffset, kernelcutoff)
 
@@ -189,10 +193,16 @@ function remorph!(self::LiftingLine, args...;
 
     calc_Dinfs!(self, Uinf)
 
+    # Reset solution
+    self.aoas .= 0
+    self.Gammas .= 0
+    self.Us .= 0
+
     return nothing
 end
 
 
+
 function save(self::LiftingLine, filename::AbstractString; 
                     format="vtk", 
                     horseshoe_suffix="_horseshoe",
@@ -265,7 +275,11 @@ function save(self::LiftingLine, filename::AbstractString;
     midpoints_data = [
                             Dict(   "field_name"  => "angleofattack",
                                     "field_type"  => "scalar",
-                                    "field_data"  => self.aoas)
+                                    "field_data"  => self.aoas),
+
+                            Dict(   "field_name"  => "U",
+                                    "field_type"  => "vector",
+                                    "field_data"  => eachcol(self.Us))
                         ]
 
     str *= gt.generateVTK(filename*midpoint_suffix, midpoints; 
@@ -288,23 +302,41 @@ function save(self::LiftingLine, filename::AbstractString;
                                 point_data=controlpoints_data, optargs...)
 
     #  ------------- OUTPUT PLANAR GEOMETRY ------------------------------------
-    planar_data = [
-                            Dict(   "field_name"  => "Gamma",
-                                    "field_type"  => "scalar",
-                                    "field_data"  => self.Gammas)
-
-                            Dict(   "field_name"  => "angleofattack",
-                                    "field_type"  => "scalar",
-                                    "field_data"  => self.aoas)
-                        ]
 
-    str *= gt.save(self.grid, filename*planar_suffix; 
-                                cell_data=planar_data, format, optargs...)
+    str *= gt.save(self.grid, filename*planar_suffix; format, optargs...)
 
     return str
 end
 
 
+function add_field(self::LiftingLine, field_name::String, field_type::String,
+                    field_data, entry_type::String;
+                    raise_warn=false, collectfield=true)
+
+    # Add field to grid
+    gt.add_field(self.grid, field_name, field_type,
+                    collectfield ? collect(field_data) : field_data, entry_type;
+                    raise_warn=raise_warn)
+
+    # Register the field
+    if !(field_name in self.fields)
+        push!(self.fields, field_name)
+    end
+
+    nothing
+end
+
+check_field(self::LiftingLine, field_name::String) = field_name in self.fields
+
+function get_field(self::LiftingLine, field_name::String)
+
+    @assert check_field(self, field_name) ""*
+        "Field $field_name not found! Available fields: $(self.fields)"
+
+    return self.grid.field[field_name]
+end
+
+
 function calc_horseshoes!(self::LiftingLine, args...; optargs...) 
     return calc_horseshoes!(self.horseshoes, 
                                 self.grid.nodes, self.linearindices, 
@@ -552,9 +584,14 @@ function solve_linear(self::LiftingLine, Uinfs::AbstractMatrix;
     # Solve system of equations
     solver(self.Gammas, self.G, self.RHS; solver_optargs...)
 
+    # Calculate velocity at lifting-line midpoints
+    self.Us .= Uinfs
+    Uind!(self, self.midpoints, self.Us)
+
     if addfields
         gt.add_field(self.grid, "Uinf", "vector", collect(eachcol(Uinfs)), "cell"; raise_warn)
         gt.add_field(self.grid, "Gamma", "scalar", self.Gammas, "cell"; raise_warn)
+        gt.add_field(self.grid, "angleofattack", "scalar", self.aoas, "cell"; raise_warn)
     end
 
     return nothing

src/FLOWPanel_liftingline_postprocess.jl

@@ -0,0 +1,166 @@
+#=##############################################################################
+# DESCRIPTION
+    Methods for postprocessing solver results of non-linear lifting line method
+
+# AUTHORSHIP
+  * Created by  : Eduardo J. Alvarez
+  * Email       : Edo.AlvarezR@gmail.com
+  * Date        : Nov 2025
+  * License     : MIT License
+=###############################################################################
+
+
+################################################################################
+# FORCE FIELDS
+################################################################################
+"""
+    calcfield_Fkj!(out::Vector, body::AbstractBody,
+                         areas::Vector, normals::Matrix, Cps::Vector,
+                         Uinf::Number, rho::Number;
+                         fieldname="F")
+
+Calculate the force of each element using the Kutta-Joukowski theorem,
+``F = \\rho \\mathbf{u} \\times \\mathbf{l} \\Gamma``. 
+``F`` is saved as a field named `fieldname`.
+
+The field is calculated in-place and added to `out` (hence, make sure that `out`
+starts with all zeroes).
+"""
+function calcfield_Fkj!(out::AbstractMatrix, ll::LiftingLine, rho::Number;
+                                                addfield=true, fieldname="Fkj")
+
+    # Error cases
+    @assert size(out, 1)==3 && size(out, 2)==ll.nelements ""*
+        "Invalid `out` matrix."*
+        " Expected size $((3, ll.nelements)); got $(size(out))."
+
+    for ei in 1:ll.nelements                # Iterate over horseshoes
+        for (i, ip1) in (                   # Iterate over bound vortices
+                            # (1, 2),           # A - Ap
+                              (2, 3),           # B - A
+                            # (3, 4)            # Bp - B
+                        )
+
+            dl1 = ll.horseshoes[1, ip1, ei] - ll.horseshoes[1, i, ei]
+            dl2 = ll.horseshoes[2, ip1, ei] - ll.horseshoes[2, i, ei]
+            dl3 = ll.horseshoes[3, ip1, ei] - ll.horseshoes[3, i, ei]
+
+            out[1, ei] += rho * ll.Gammas[ei] * (ll.Us[2, ei]*dl3 - ll.Us[3, ei]*dl2)
+            out[2, ei] += rho * ll.Gammas[ei] * (ll.Us[3, ei]*dl1 - ll.Us[1, ei]*dl3)
+            out[3, ei] += rho * ll.Gammas[ei] * (ll.Us[1, ei]*dl2 - ll.Us[2, ei]*dl1)
+
+        end
+    end
+
+    # Save field in lifting line
+    if addfield
+        add_field(ll, fieldname, "vector", eachcol(out), "cell")
+    end
+
+    return out
+end
+
+
+"""
+    calcfield_Fkj(args...; optargs...)
+
+Similar to [`calcfield_Fkj!`](@ref) but automatically pre-allocating `out` if it
+hasn't been pre-allocated yet
+"""
+function calcfield_Fkj(ll::LiftingLine{R}, args...; fieldname="Fkj", optargs...) where {R}
+    
+    if fieldname in ll.fields
+        out = get_field(ll, fieldname)["field_data"]
+        out .= 0
+    else
+        out = zeros(R, 3, ll.nelements)
+    end
+
+    return calcfield_Fkj!(out, ll, args...; optargs...)
+
+end
+
+
+
+
+"""
+    calcfield_fkj!(out::Vector, body::AbstractBody,
+                         areas::Vector, normals::Matrix, Cps::Vector,
+                         Uinf::Number, rho::Number;
+                         fieldname="F")
+
+Calculate loading distribution (force per unit span) using the Kutta-Joukowski 
+theorem,
+``f = \\rho \\mathbf{u} \\times \\hat{\\mathbf{y}} \\Gamma``. 
+``f`` is saved as a field named `fieldname`.
+
+The field is calculated in-place and added to `out` (hence, make sure that `out`
+starts with all zeroes).
+"""
+function calcfield_fkj!(out::AbstractMatrix, ll::LiftingLine, rho::Number;
+                                                addfield=true, fieldname="fkj")
+
+    # Error cases
+    @assert size(out, 1)==3 && size(out, 2)==ll.nelements ""*
+        "Invalid `out` matrix."*
+        " Expected size $((3, ll.nelements)); got $(size(out))."
+
+    for ei in 1:ll.nelements                # Iterate over horseshoes
+        for (i, ip1) in (                   # Iterate over bound vortices
+                            # (1, 2),           # A - Ap
+                              (2, 3),           # B - A
+                            # (3, 4)            # Bp - B
+                        )
+
+            # Filament length
+            dl1 = ll.horseshoes[1, ip1, ei] - ll.horseshoes[1, i, ei]
+            dl2 = ll.horseshoes[2, ip1, ei] - ll.horseshoes[2, i, ei]
+            dl3 = ll.horseshoes[3, ip1, ei] - ll.horseshoes[3, i, ei]
+
+            # Velocity
+            U1 = ll.Us[1, ei]
+            U2 = ll.Us[2, ei]
+            U3 = ll.Us[3, ei]
+            magU = sqrt(U1^2 + U2^2 + U3^2)
+
+            # Calculate tranversal length (counter-projection of U): dy = dl - (dl⋅hatU)hatU
+            dldotU = dl1*U1 + dl2*U2 + dl3*U3
+            dy1 = dl1 - dldotU*U1/magU^2
+            dy2 = dl2 - dldotU*U2/magU^2
+            dy3 = dl3 - dldotU*U3/magU^2
+            magdy = sqrt(dy1^2 + dy2^2 + dy3^2)
+
+            out[1, ei] += rho * ll.Gammas[ei] * (U2*dy3 - U3*dy2)/magdy
+            out[2, ei] += rho * ll.Gammas[ei] * (U3*dy1 - U1*dy3)/magdy
+            out[3, ei] += rho * ll.Gammas[ei] * (U1*dy2 - U2*dy1)/magdy
+
+        end
+    end
+
+    # Save field in lifting line
+    if addfield
+        add_field(ll, fieldname, "vector", eachcol(out), "cell")
+    end
+
+    return out
+end
+
+
+"""
+    calcfield_fkj(args...; optargs...)
+
+Similar to [`calcfield_fkj!`](@ref) but automatically pre-allocating `out` if it
+hasn't been pre-allocated yet
+"""
+function calcfield_fkj(ll::LiftingLine{R}, args...; fieldname="fkj", optargs...) where {R}
+    
+    if fieldname in ll.fields
+        out = get_field(ll, fieldname)["field_data"]
+        out .= 0
+    else
+        out = zeros(R, 3, ll.nelements)
+    end
+
+    return calcfield_fkj!(out, ll, args...; optargs...)
+
+end

src/FLOWPanel_liftingline_stripwise.jl

@@ -13,7 +13,7 @@ import CSV
 import DataFrames: DataFrame
 import AirfoilPrep
 
-abstract type StripwiseElement end
+abstract type StripwiseElement <: AbstractElement end
 
 ################################################################################
 # JETFOIL ELEMENT STRUCT