Add a compact-support core regularization function for OLAF particle methods

docs/source/user/aerodyn-olaf/ExampleFiles/ExampleFile--OLAF.dat

@@ -22,7 +22,8 @@ default FWShedVorticity    - Include shed vorticity in the far wake {default: fa
 ------------------- WAKE REGULARIZATIONS AND DIFFUSION -----------------------------------------
 default DiffusionMethod    - Diffusion method to account for viscous effects {0: None, 1: Core Spreading, "default": 0}         
 2       RegDeterMethod     - Method to determine the regularization parameters {0: Constant, 1: Optimized, 2: Chord-scaled, 3: dr-scaled, default: 0 }          
-default RegFunction        - Viscous diffusion function {0: None, 1: Rankine, 2: LambOseen, 3: Vatistas, 4: Denominator, "default": 3} (switch)
+default RegFunction        - Segment regularization function {0: None, 1: Rankine, 2: LambOseen, 3: Vatistas, 4: Denominator, "default": 3} (switch)
+default RegFunctionPart    - Particle regularization function {0: None, 1: Exponential, 2: Compact, "default": 1} [only if VelocityMethod=2,3] (switch)
 default WakeRegMethod      - Wake regularization method {1: Constant, 2: Stretching, 3: Age, default: 1} (switch)
 0.25    WakeRegFactor      - Wake regularization factor (m or -)
 0.25    WingRegFactor      - Wing regularization factor (m or -)

docs/source/user/aerodyn-olaf/InputFiles.rst

@@ -178,13 +178,29 @@ See :numref:`Guidelines-OLAF` for recommendations on setting up this parameter.
 
 
 **RegFunction** [switch] specifies the regularization function used to remove
-the singularity of the vortex elements, as specified in
+the singularity of the vortex elements/segments, as specified in
 :numref:`sec:vortconv`. There are five options: 1) no correction *[0]*,
 2) the Rankine method *[1]*, 3) the Lamb-Oseen method *[2]*, 4) the Vatistas
 method *[3]*, and 5) the denominator offset method *[4]*. 
 The functions are given in :numref:`sec:RegularizationFunction`. 
 The default option is *[3]*.
 
+**RegFunctionPart** [switch] specifies the regularization function used for the
+vortex particles, which are used when a particle-based velocity method is
+selected (*VelocityMethod* = *[2,3]*). There are three options: 1) no
+correction *[0]*, 2) the exponential method *[1]*, and 3) the compact-support
+method *[2]*. The functions are given in
+:numref:`sec:RegularizationFunctionPart`. 
+The compact-support option *[2]* has a finite support radius beyond which the
+kernel reverts to the exact singular kernel; with the tree-accelerated particle
+method (*VelocityMethod* = *[2]*) this allows more particle interactions to be
+handled by the far-field multipole approximation, so it can be significantly
+faster than the exponential option depending on the case (wake size, particle
+count, and regularization parameter).
+The compact-support kernel is also purely polynomial and avoids evaluating a
+transcendental (exponential) function, which can further reduce cost.
+The default option is *[1]*.
+
 **WakeRegMethod** [switch] specifies the method of determining viscous core
 radius (i.e., the regularization parameter). There are three options: 1)
 constant *[1]*, 2) stretching *[2]*, and 3) age *[3]*. The methods are

docs/source/user/aerodyn-olaf/OLAFTheory.rst

@@ -489,8 +489,16 @@ refinement of this option will be considered in the future.
 
 .. _sec:RegularizationFunction:
 
-Implemented regularization functions
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+Segment regularization functions
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The regularization functions described in this section apply to the vortex
+*segments* used to represent both the bound (blade) vorticity and the wake
+vorticity, i.e., they regularize the segment Biot-Savart kernel of
+Eq. :eq:`eq:BiotSavartSegment`. They are selected with the input
+**RegFunction**. The corresponding regularization functions used by the
+vortex-particle representation of the wake are described separately in
+:numref:`sec:RegularizationFunctionPart`.
 
 Several regularization functions have been
 developed (:cite:`olaf-Rankine58_1,olaf-Scully75_1,olaf-Vatistas91_1`).  At present, five
@@ -563,6 +571,109 @@ denominator of Eq. :eq:`eq:BiotSavartSegment`, proportional to the filament
 length :math:`r_0`. In this case, :math:`F_\nu=1`. This method is found in the
 work of van Garrel (:cite:`olaf-Garrel03_1`).
 
+.. _sec:RegularizationFunctionPart:
+
+Particle regularization functions
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The regularization functions of :numref:`sec:RegularizationFunction` apply to
+the vortex segments. When a particle-based velocity method is selected
+(**VelocityMethod=[3]**, the direct vortex-particle method, or
+**VelocityMethod=[2]**, the tree-accelerated particle method), the wake segments
+are converted into vortex particles and a separate set of regularization
+functions applies. The particle regularization function is selected with the
+input **RegFunctionPart**.
+
+Segment-to-particle conversion
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Each wake segment of circulation :math:`\Gamma` and length vector
+:math:`\vec{l}=\vec{x}_2-\vec{x}_1` is divided into :math:`n_p`
+(**PartPerSegment**) equally spaced particles. The intensity of a particle
+:math:`\vec{\alpha}` (the vorticity integrated over the volume represented by
+the particle, :math:`\vec{\alpha}=\vec{\omega}\,dV`) is obtained from the parent
+segment as
+
+.. math::
+   \vec{\alpha} = \frac{\Gamma\,\vec{l}}{n_p}
+   :label: eq:SegToPart
+
+so that the :math:`n_p` particles of a segment sum to the total segment
+vorticity :math:`\Gamma\,\vec{l}`. The particles are placed at the centers of
+the :math:`n_p` sub-segments, and each particle inherits the regularization
+parameter (core radius :math:`r_c`) of its parent segment.
+
+Particle velocity kernel
+^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The velocity induced at a point :math:`\vec{x}` by a vortex particle of
+intensity :math:`\vec{\alpha}` located at :math:`\vec{x}_p` is the regularized
+point-vortex kernel
+
+.. math::
+   \vec{v}(\vec{x}) = \frac{1}{4\pi}\, g(\bar{r})\,
+   \frac{\vec{\alpha}\times\vec{r}}{r^3}
+   ,\qquad \vec{r}=\vec{x}-\vec{x}_p,\quad r=|\vec{r}|,\quad \bar{r}=\frac{r}{r_c}
+   :label: eq:BiotSavartParticle
+
+where :math:`g` is the particle regularization factor (the particle-method
+analog of :math:`F_\nu`) and :math:`r_c` is the particle core radius. Away from
+the particle center the induced velocity decays as :math:`1/r^2`. If no
+correction is used (**RegFunctionPart=[0]**), :math:`g=1` and the singular
+point-vortex kernel is recovered.
+
+Exponential
+^^^^^^^^^^^
+
+If the exponential method is used (**RegFunctionPart=[1]**, the default), the
+regularization factor is
+
+.. math::
+   g(\bar{r}) = 1 - \exp\!\left(-\bar{r}^{\,3}\right)
+   :label: eq:PartExp
+
+For :math:`\bar{r} > 2` (beyond two core radii) the exponential term is
+negligible and :math:`g` is set to :math:`1`.
+
+Compact support
+^^^^^^^^^^^^^^^^
+
+If the compact-support method is used (**RegFunctionPart=[2]**), the kernel is
+*exactly* singular (:math:`g=1`) beyond a finite support radius
+:math:`r_c'=1.6\,r_c`, and inside the support it is given by a polynomial
+mollifier
+
+.. math::
+   g(\bar{r}') = \begin{cases}
+   \dfrac{35\,\bar{r}'^{3} - 42\,\bar{r}'^{5} + 15\,\bar{r}'^{7}}{8}
+   & 0 \le \bar{r}' < 1 \\[2mm]
+   1 & \bar{r}' \ge 1
+   \end{cases}
+   ,\qquad \bar{r}' = \frac{r}{r_c'} = \frac{r}{1.6\,r_c}
+   :label: eq:PartCompact
+
+The polynomial is smooth (:math:`C^\infty`) inside the support; the overall
+regularization factor is globally :math:`C^2`, its regularity being limited by
+the match to the singular kernel at the support boundary, where :math:`g=1`,
+:math:`g'=0`, and :math:`g''=0` but :math:`g'''\neq 0` at :math:`\bar{r}'=1`. At
+the center the mollifier vanishes as :math:`\bar{r}'^{3}`
+(:math:`g=g'=g''=0`), so the induced velocity is regular there. The
+support-radius factor :math:`1.6` makes the compact kernel roughly equivalent
+to the exponential kernel with the same :math:`r_c`. It matches the peak
+tangential (swirl) velocity *magnitude* of the two kernels, while being a
+reasonable compromise across three possible equivalence criteria:
+
+.. math::
+   \frac{r_c'}{r_c} \approx \begin{cases}
+   1.54 & \text{radius of peak tangential velocity} \\
+   1.60 & \text{peak tangential velocity magnitude} \\
+   1.65 & \text{second moment of the vorticity distribution}
+   \end{cases}
+
+Because the kernel is exactly equal to the singular kernel beyond :math:`r_c'`,
+the tree-accelerated particle method can apply an exact far-field cutoff at
+that radius.
+
 .. _sec:corerad:
 
 Time Evolution of the Regularization Parameter–Core Spreading Method

docs/source/user/api_change.rst

@@ -16,6 +16,8 @@ Under-relaxation is introduced for the tight-coupling iterative solver to improv
 
 The generalized support-structure (GS) influence model was added to AeroDyn. This introduces three new switches (``GSPotent``, ``GSShadow``, ``GSAero``) after ``TwrAero`` in the AeroDyn primary input file, and two new sections (``General support structure joints`` and ``General support structure members``) after the ``Tower Influence and Aerodynamics`` section. The two new sections are required even when the GS model is disabled (set ``NumGSJoints`` and ``NumGSMembers`` to 0, keeping the table header lines).
 
+OLAF now selects the regularization function for the vortex particles separately from the vortex segments. The new ``RegFunctionPart`` input on line 26 of the OLAF input file controls the particle kernel, and ``RegFunction`` now applies only to the segments. Previously the particle kernel was inferred from ``RegFunction``: any regularized value gave an exponential particle kernel, while ``RegFunction=0`` gave an unregularized one. The default ``RegFunctionPart=1`` reproduces the exponential kernel, so results are unchanged for decks with ``RegFunction`` greater than 0. Decks with ``RegFunction=0`` will change: the particles are now regularized with the exponential kernel unless ``RegFunctionPart=0`` is also set. Setting ``RegFunctionPart=0`` restores an unregularized particle kernel, but with ``VelocityMethod=2`` the results will still differ from previous versions because the particle-tree far-field cutoff is no longer padded by the particle core radius when the particle kernel is unregularized.
+
 ============================================= ======== ==================== ==========================================================================================================================================================================================================================================================================================================
 Added in OpenFAST `5.1.0`
 --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
@@ -36,6 +38,7 @@ AeroDyn                                       \*                            ====
 AeroDyn                                       \*       NumGSMembers         0   NumGSMembers - Number of general support members (-)
 AeroDyn                                       \*                            GSMemberID  GSMJointID1  GSMJointID2     GSMDia1       GSMDia2     GSMCd1        GSMCd2        GSMTI1       GSMTI2       GSMDiv
 AeroDyn                                       \*                            (-)            (-)          (-)           (m)           (m)          (-)           (-)          (-)          (-)           (m)
+OLAF                                          26       RegFunctionPart            2  RegFunctionPart - Particle regularization function {0: None, 1: Exponential, 2: Compact, "default": 1} [only if VelocityMethod=2,3] (switch)
 ============================================= ======== ==================== ==========================================================================================================================================================================================================================================================================================================
 
 OpenFAST v4.2.x to OpenFAST  v5.0.0

modules/aerodyn/src/FVW.f90

@@ -453,6 +453,7 @@ SUBROUTINE FVW_SetParametersFromInputFile( InputFileData, p, ErrStat, ErrMsg )
    p%FWShedVorticity      = InputFileData%FWShedVorticity
    p%DiffusionMethod      = InputFileData%DiffusionMethod
    p%RegFunction          = InputFileData%RegFunction
+   p%RegFunctionPart      = InputFileData%RegFunctionPart
    p%RegDeterMethod       = InputFileData%RegDeterMethod
    p%WakeRegMethod        = InputFileData%WakeRegMethod
    p%WakeRegParam         = InputFileData%WakeRegParam

modules/aerodyn/src/FVW_BiotSavart.f90

@@ -12,6 +12,7 @@ module FVW_BiotSavart
    real(ReKi),parameter :: MIN_EXP_VALUE=-10.0_ReKi
    real(ReKi),parameter :: PART_REG_NRAD =  2.0_ReKi             !< Particle exp mollifier treated as 1 beyond this many core radii (matches 2*rc far-field multipole floor)
    real(ReKi),parameter :: PART_REG_CUT3 =  PART_REG_NRAD**3     !< Corresponding (r/rc)^3 cutoff
+   real(ReKi),parameter :: PART_REG_C2   =  1.6_ReKi             !< Compact-C2 support = PART_REG_C2*RegParam (exp-equivalent core); also the exact multipole floor for that kernel
    real(ReKi),parameter :: MINDENOM=0.0_ReKi
 !    real(ReKi),parameter :: MINDENOM=1e-15_ReKi
    real(ReKi),parameter :: MINNORM=1e-4
@@ -31,6 +32,21 @@ module FVW_BiotSavart
 
 contains
 
+!> Multipole floor factor for the particle kernels: number of core radii beyond which the
+!! regularized kernel equals the singular 1/r^3, so the far-field multipole expansion is valid.
+pure function PartRegFloorFactor(RegFunction) result(f)
+   integer(IntKi), intent(in) :: RegFunction
+   real(ReKi)                 :: f
+   select case (RegFunction)
+   case (idRegNone)    ! Unregularized kernel: no extra core-based floor needed
+      f = 0.0_ReKi
+   case (idRegCompact) ! Truly compact: exactly singular beyond rc = PART_REG_C2*RegParam
+      f = PART_REG_C2
+   case default        ! Exponential: conservative 2*rc floor
+      f = PART_REG_NRAD
+   end select
+end function PartRegFloorFactor
+
 
 !> Induced velocity from one segment at one control points
 subroutine ui_seg_11(DeltaPa, DeltaPb, SegGamma, RegFunction, RegParam1, Uind)
@@ -376,7 +392,8 @@ subroutine ui_part_nograd_11(DeltaP, Alpha, RegFunction, RegParam, Ui)
    real(ReKi),dimension(3) :: C          !< Cross product of Alpha and r
    real(ReKi)              :: E          !< Exponential poart for the mollifider
    real(ReKi)              :: r3_inv     !< 
-   real(ReKi)              :: r2, r3, rc3!< |r|^2, |r|^3, RegParam^3 (reused to avoid recomputing ** intrinsics)
+   real(ReKi)              :: r2, r3, rc3!< |r|^2, |r|^3, core^3 (reused to avoid recomputing ** intrinsics)
+   real(ReKi)              :: rc2, t     !< compact-C2 core^2 and (r/rc)^2
    real(ReKi)              :: rDeltaP    !< norm , distance between point and particle
    real(ReKi)              :: ScalarPart !< the part containing the inverse of the distance, but not 4pi, Mollifier
    r2      = DeltaP(1)**2+ DeltaP(2)**2+ DeltaP(3)**2
@@ -402,10 +419,16 @@ subroutine ui_part_nograd_11(DeltaP, Alpha, RegFunction, RegParam, Ui)
             E          = exp(-r3/rc3)
             ScalarPart = (1._ReKi-E)*r3_inv*fourpi_inv
          endif
-      case (idRegCompact) ! Compact support
-         rc3        = RegParam*RegParam*RegParam
-         r3_inv     = 1._ReKi/sqrt(rc3*rc3+r3*r3)
-         ScalarPart = r3_inv*fourpi_inv
+      case (idRegCompact) ! Truly compact C2 core: exactly singular for r>=rc, rc=PART_REG_C2*RegParam
+         rc2 = (PART_REG_C2*RegParam)**2
+         if (r2 >= rc2) then ! outside core: kernel is exactly the singular 1/r^3
+            r3_inv     = 1._ReKi/r3
+            ScalarPart = r3_inv*fourpi_inv
+         else                ! inside core: C2 polynomial mollifier
+            t          = r2/rc2
+            rc3        = rc2*PART_REG_C2*RegParam
+            ScalarPart = (35._ReKi + t*(-42._ReKi + 15._ReKi*t))*fourpi_inv/(8._ReKi*rc3)
+         endif
       case default 
          print*,'[ERROR] Wrong regularization function for particles',RegFunction
          STOP

modules/aerodyn/src/FVW_IO.f90

@@ -60,6 +60,7 @@ SUBROUTINE FVW_ReadInputFile( FileName, p, m, Inp, ErrStat, ErrMsg )
    CALL ReadVarWDefault(UnIn,FileName,Inp%DiffusionMethod    ,'DiffusionMethod'   ,'',idDiffusionNone , ErrStat2,ErrMsg2); if(Failed())return
    CALL ReadVarWDefault(UnIn,FileName,Inp%RegDeterMethod     ,'RegDeterMethod'    ,'',idRegDeterConstant, ErrStat2,ErrMsg2); if(Failed())return
    CALL ReadVarWDefault(UnIn,FileName,Inp%RegFunction        ,'RegFunction'       ,'',idRegVatistas   , ErrStat2,ErrMsg2); if(Failed())return
+   CALL ReadVarWDefault(UnIn,FileName,Inp%RegFunctionPart    ,'RegFunctionPart'   ,'',idRegExp        , ErrStat2,ErrMsg2); if(Failed())return
    CALL ReadVarWDefault(UnIn,FileName,Inp%WakeRegMethod      ,'WakeRegMethod'     ,'',idRegAge        , ErrStat2,ErrMsg2); if(Failed())return
    CALL ReadVar        (UnIn,FileName,Inp%WakeRegParam       ,'WakeRegParam'      ,''                 , ErrStat2,ErrMsg2); if(Failed())return
    CALL ReadVar        (UnIn,FileName,Inp%WingRegParam       ,'WingRegParam'      ,''                 , ErrStat2,ErrMsg2); if(Failed())return
@@ -191,6 +192,7 @@ SUBROUTINE FVW_ReadInputFile( FileName, p, m, Inp, ErrStat, ErrMsg )
    if (Check(.not.(ANY(idDiffusionVALID==Inp%DiffusionMethod)) , 'Diffusion method (DiffusionMethod) not implemented: '//trim(Num2LStr(Inp%DiffusionMethod)))) return
    if (Check(.not.(ANY(idRegDeterVALID ==Inp%RegDeterMethod))  , 'Regularization determination method (RegDeterMethod) not yet implemented: '//trim(Num2LStr(Inp%RegDeterMethod)))) return
    if (Check(.not.(ANY(idRegVALID      ==Inp%RegFunction  )), 'Regularization function (RegFunction) not implemented: '//trim(Num2LStr(Inp%RegFunction)))) return
+   if (Check(.not.(ANY(idRegPartVALID  ==Inp%RegFunctionPart)), 'Particle regularization function (RegFunctionPart) not implemented: '//trim(Num2LStr(Inp%RegFunctionPart)))) return
    if (Check(.not.(ANY(idRegMethodVALID==Inp%WakeRegMethod)), 'Wake regularization method (WakeRegMethod) not implemented: '//trim(Num2LStr(Inp%WakeRegMethod)))) return
    if (Check(.not.(ANY(idShearVALID    ==Inp%ShearModel   )), 'Shear model (ShearModel) not valid: '//trim(Num2LStr(Inp%ShearModel)))) return
    if (Check(.not.(ANY(idVelocityVALID ==Inp%VelocityMethod(1))), 'Velocity method (VelocityMethod(1)) not valid: '//trim(Num2LStr(Inp%VelocityMethod(1))))) return