diff --git a/docs/source/user/aerodyn-olaf/ExampleFiles/ExampleFile--OLAF.dat b/docs/source/user/aerodyn-olaf/ExampleFiles/ExampleFile--OLAF.dat index 4fdfb0ad97..2763378554 100644 --- a/docs/source/user/aerodyn-olaf/ExampleFiles/ExampleFile--OLAF.dat +++ b/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 -) diff --git a/docs/source/user/aerodyn-olaf/InputFiles.rst b/docs/source/user/aerodyn-olaf/InputFiles.rst index afc81f7b85..19c97b188e 100644 --- a/docs/source/user/aerodyn-olaf/InputFiles.rst +++ b/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 diff --git a/docs/source/user/aerodyn-olaf/OLAFTheory.rst b/docs/source/user/aerodyn-olaf/OLAFTheory.rst index ac21622da4..1379f4e626 100644 --- a/docs/source/user/aerodyn-olaf/OLAFTheory.rst +++ b/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 diff --git a/docs/source/user/api_change.rst b/docs/source/user/api_change.rst index f877be8949..97040f95bd 100644 --- a/docs/source/user/api_change.rst +++ b/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 diff --git a/modules/aerodyn/src/FVW.f90 b/modules/aerodyn/src/FVW.f90 index 15da86e339..ba78825ab9 100644 --- a/modules/aerodyn/src/FVW.f90 +++ b/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 diff --git a/modules/aerodyn/src/FVW_BiotSavart.f90 b/modules/aerodyn/src/FVW_BiotSavart.f90 index d6db453b6f..bf93782203 100644 --- a/modules/aerodyn/src/FVW_BiotSavart.f90 +++ b/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 diff --git a/modules/aerodyn/src/FVW_IO.f90 b/modules/aerodyn/src/FVW_IO.f90 index aef35dcfcb..b8bfb615d6 100644 --- a/modules/aerodyn/src/FVW_IO.f90 +++ b/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 diff --git a/modules/aerodyn/src/FVW_Registry.txt b/modules/aerodyn/src/FVW_Registry.txt index 2d9ce90812..72968d627f 100644 --- a/modules/aerodyn/src/FVW_Registry.txt +++ b/modules/aerodyn/src/FVW_Registry.txt @@ -116,6 +116,7 @@ typedef ^ ^ IntKi typedef ^ ^ ReKi CoreSpreadEddyVisc - - - "Eddy viscosity used in the core spreading method" typedef ^ ^ IntKi RegDeterMethod - - - "Regularization determinatino method (manual, automatic)" - typedef ^ ^ IntKi RegFunction - - - "Type of regularizaion function (LambOseen, Vatistas, see FVW_BiotSavart)" - +typedef ^ ^ IntKi RegFunctionPart - - - "Type of particle regularization function (None, Exp, Compact; see FVW_BiotSavart)" - typedef ^ ^ IntKi WakeRegMethod - - - "Method for regularization (constant, stretching, age, etc.)" - typedef ^ ^ ReKi WakeRegParam - - - "Initial value of the regularization parameter" typedef ^ ^ ReKi WingRegParam - - - "Regularization parameter of the wing" @@ -334,6 +335,7 @@ typedef ^ ^ IntKi typedef ^ ^ ReKi CoreSpreadEddyVisc - - - "Eddy viscosity used in the core spreading method" typedef ^ ^ IntKi RegDeterMethod - - - "Regularization determinatino method (manual, automatic)" - typedef ^ ^ IntKi RegFunction - - - "Type of regularizaion function (LambOseen, Vatistas, see FVW_BiotSavart)" - +typedef ^ ^ IntKi RegFunctionPart - - - "Type of particle regularization function (None, Exp, Compact; see FVW_BiotSavart)" - typedef ^ ^ IntKi WakeRegMethod - - - "Method for regularization (constant, stretching, age, etc.)" - typedef ^ ^ ReKi WakeRegParam - - - "Factor used in the regularization " typedef ^ ^ ReKi WingRegParam - - - "Factor used in the regularization " diff --git a/modules/aerodyn/src/FVW_Subs.f90 b/modules/aerodyn/src/FVW_Subs.f90 index 43ae7d3eb5..95c3aa0b2c 100644 --- a/modules/aerodyn/src/FVW_Subs.f90 +++ b/modules/aerodyn/src/FVW_Subs.f90 @@ -843,7 +843,7 @@ subroutine FVW_InitMiscVarsPostParam( p, m, ErrStat, ErrMsg ) call AllocAry( m%Part%Alpha , 3, nPart, 'PartAlpha' , ErrStat2, ErrMsg2 ); if(Failed())return; m%Part%Alpha = -999999_ReKi; call AllocAry( m%Part%RegParam, nPart, 'PartEpsilon', ErrStat2, ErrMsg2 ); if(Failed())return; m%Part%RegParam= -999999_ReKi; m%Part%nAct = -1 ! Active particles - m%Part%RegFunction = p%RegFunction + m%Part%RegFunction = p%RegFunctionPart endif ! TODO Figure out Uind, CPs needed for grid @@ -1185,11 +1185,11 @@ subroutine InducedVelocitiesAll_OnGrid_Calc(g, p, Sgmt, Part, Tree, Panl, ErrSta end subroutine InducedVelocitiesAll_OnGrid_Calc !> Wrapper to setup part from set of segments -subroutine SegmentsToPartWrap(Sgmt, nSeg, PartPerSegment, RegFunction, Part, allocPart) +subroutine SegmentsToPartWrap(Sgmt, nSeg, PartPerSegment, RegFunctionPart, Part, allocPart) type(T_Sgmt), intent(in ) :: Sgmt !< Segments integer(IntKi), intent(in ) :: nSeg !< Number of segments to use (might not use all of them) integer(IntKi), intent(in ) :: PartPerSegment !< Number of particles per segment - integer(IntKi), intent(in ) :: RegFunction !< Regularization function + integer(IntKi), intent(in ) :: RegFunctionPart !< Particle regularization function (None, Exp, Compact) type(T_Part), intent(inout) :: Part !< Particles logical, intent(in ) :: allocPart !< allocate particles integer(IntKi) :: iHeadP @@ -1220,9 +1220,7 @@ subroutine SegmentsToPartWrap(Sgmt, nSeg, PartPerSegment, RegFunction, Part, all Part%nAct = nPart ! TODO add iHeadPart if particles already present call SegmentsToPart(Sgmt%Points, Sgmt%Connct, Sgmt%Gamma, Sgmt%Epsilon, 1, nSeg, PartPerSegment, Part%P, Part%Alpha, Part%RegParam, iHeadP) - if (RegFunction/=idRegNone) then - Part%RegFunction = idRegExp ! TODO need to find a good equivalence and potentially adapt Epsilon in SegmentsToPart - endif + Part%RegFunction = RegFunctionPart ! TODO Epsilon equivalence between segment and particle cores may still need tuning in SegmentsToPart if (DEV_VERSION) then call find_nan_2D(Part%P(:,1:nPart) , 'SegmentsToPartWrap Part%P') call find_nan_2D(Part%Alpha(:,1:nPart), 'SegmentsToPartWrap Part%Alpha') @@ -1265,10 +1263,10 @@ subroutine InducedVelocitiesAll_Init(p, x, m, Sgmt, Part, Tree, Panl, ErrStat, E select case (p%VelocityMethod(iVel)) case (idVelocityPart) ! --- Convert segments to particles - call SegmentsToPartWrap(Sgmt, nSeg, p%PartPerSegment(iVel), p%RegFunction, Part, allocPart=allocPart) + call SegmentsToPartWrap(Sgmt, nSeg, p%PartPerSegment(iVel), p%RegFunctionPart, Part, allocPart=allocPart) case (idVelocityTreePart) ! --- Convert segments to particles - call SegmentsToPartWrap(Sgmt, nSeg, p%PartPerSegment(iVel), p%RegFunction, Part, allocPart=allocPart) + call SegmentsToPartWrap(Sgmt, nSeg, p%PartPerSegment(iVel), p%RegFunctionPart, Part, allocPart=allocPart) ! --- Grow particle tree call grow_tree_part(Tree, Part%nAct, Part%P, Part%Alpha, Part%RegFunction, Part%RegParam, 0) case (idVelocityTreeSeg) @@ -1566,7 +1564,7 @@ subroutine LiftingLineInducedVelocities(p, x, InductionAtCP, iDepthStart, m, Err call ui_seg( 1, nCPs, CPs, 1, nSeg, m%Sgmt%Points, m%Sgmt%Connct, m%Sgmt%Gamma, m%Sgmt%RegFunction, m%Sgmt%Epsilon, Uind) case (idVelocityPart) - call SegmentsToPartWrap(m%Sgmt, nSeg, p%PartPerSegment(iVel), p%RegFunction, m%Part, allocPart=.false.) + call SegmentsToPartWrap(m%Sgmt, nSeg, p%PartPerSegment(iVel), p%RegFunctionPart, m%Part, allocPart=.false.) call ui_part_nograd(nCPs, CPs, m%Part%nAct, m%Part%P, m%Part%Alpha, m%Part%RegFunction, m%Part%RegParam, Uind) !deallocate(Part%P, Part%Alpha, Part%RegParam) @@ -1576,7 +1574,7 @@ subroutine LiftingLineInducedVelocities(p, x, InductionAtCP, iDepthStart, m, Err call cut_tree(Tree) case (idVelocityTreePart) - call SegmentsToPartWrap(m%Sgmt, nSeg, p%PartPerSegment(iVel), p%RegFunction, m%Part, allocPart=.false.) + call SegmentsToPartWrap(m%Sgmt, nSeg, p%PartPerSegment(iVel), p%RegFunctionPart, m%Part, allocPart=.false.) call grow_tree_part(Tree, m%Part%nAct, m%Part%P, m%Part%Alpha, m%Part%RegFunction, m%Part%RegParam, 0) call ui_tree_part(Tree, nCPs, CPs, p%TreeBranchFactor(iVel), DistanceDirect, Uind, ErrStat, ErrMsg) !deallocate(Part%P, Part%Alpha, Part%RegParam) diff --git a/modules/aerodyn/src/FVW_Types.f90 b/modules/aerodyn/src/FVW_Types.f90 index 9f0b0b1aa7..e8d100f37b 100644 --- a/modules/aerodyn/src/FVW_Types.f90 +++ b/modules/aerodyn/src/FVW_Types.f90 @@ -157,6 +157,7 @@ MODULE FVW_Types REAL(ReKi) :: CoreSpreadEddyVisc = 0.0_ReKi !< Eddy viscosity used in the core spreading method [-] INTEGER(IntKi) :: RegDeterMethod = 0_IntKi !< Regularization determinatino method (manual, automatic) [-] INTEGER(IntKi) :: RegFunction = 0_IntKi !< Type of regularizaion function (LambOseen, Vatistas, see FVW_BiotSavart) [-] + INTEGER(IntKi) :: RegFunctionPart = 0_IntKi !< Type of particle regularization function (None, Exp, Compact; see FVW_BiotSavart) [-] INTEGER(IntKi) :: WakeRegMethod = 0_IntKi !< Method for regularization (constant, stretching, age, etc.) [-] REAL(ReKi) :: WakeRegParam = 0.0_ReKi !< Initial value of the regularization parameter [-] REAL(ReKi) :: WingRegParam = 0.0_ReKi !< Regularization parameter of the wing [-] @@ -385,6 +386,7 @@ MODULE FVW_Types REAL(ReKi) :: CoreSpreadEddyVisc = 0.0_ReKi !< Eddy viscosity used in the core spreading method [-] INTEGER(IntKi) :: RegDeterMethod = 0_IntKi !< Regularization determinatino method (manual, automatic) [-] INTEGER(IntKi) :: RegFunction = 0_IntKi !< Type of regularizaion function (LambOseen, Vatistas, see FVW_BiotSavart) [-] + INTEGER(IntKi) :: RegFunctionPart = 0_IntKi !< Type of particle regularization function (None, Exp, Compact; see FVW_BiotSavart) [-] INTEGER(IntKi) :: WakeRegMethod = 0_IntKi !< Method for regularization (constant, stretching, age, etc.) [-] REAL(ReKi) :: WakeRegParam = 0.0_ReKi !< Factor used in the regularization [-] REAL(ReKi) :: WingRegParam = 0.0_ReKi !< Factor used in the regularization [-] @@ -1541,6 +1543,7 @@ subroutine FVW_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg) DstParamData%CoreSpreadEddyVisc = SrcParamData%CoreSpreadEddyVisc DstParamData%RegDeterMethod = SrcParamData%RegDeterMethod DstParamData%RegFunction = SrcParamData%RegFunction + DstParamData%RegFunctionPart = SrcParamData%RegFunctionPart DstParamData%WakeRegMethod = SrcParamData%WakeRegMethod DstParamData%WakeRegParam = SrcParamData%WakeRegParam DstParamData%WingRegParam = SrcParamData%WingRegParam @@ -1641,6 +1644,7 @@ subroutine FVW_PackParam(RF, Indata) call RegPack(RF, InData%CoreSpreadEddyVisc) call RegPack(RF, InData%RegDeterMethod) call RegPack(RF, InData%RegFunction) + call RegPack(RF, InData%RegFunctionPart) call RegPack(RF, InData%WakeRegMethod) call RegPack(RF, InData%WakeRegParam) call RegPack(RF, InData%WingRegParam) @@ -1719,6 +1723,7 @@ subroutine FVW_UnPackParam(RF, OutData) call RegUnpack(RF, OutData%CoreSpreadEddyVisc); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%RegDeterMethod); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%RegFunction); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%RegFunctionPart); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%WakeRegMethod); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%WakeRegParam); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%WingRegParam); if (RegCheckErr(RF, RoutineName)) return @@ -4232,6 +4237,7 @@ subroutine FVW_CopyInputFile(SrcInputFileData, DstInputFileData, CtrlCode, ErrSt DstInputFileData%CoreSpreadEddyVisc = SrcInputFileData%CoreSpreadEddyVisc DstInputFileData%RegDeterMethod = SrcInputFileData%RegDeterMethod DstInputFileData%RegFunction = SrcInputFileData%RegFunction + DstInputFileData%RegFunctionPart = SrcInputFileData%RegFunctionPart DstInputFileData%WakeRegMethod = SrcInputFileData%WakeRegMethod DstInputFileData%WakeRegParam = SrcInputFileData%WakeRegParam DstInputFileData%WingRegParam = SrcInputFileData%WingRegParam @@ -4281,6 +4287,7 @@ subroutine FVW_PackInputFile(RF, Indata) call RegPack(RF, InData%CoreSpreadEddyVisc) call RegPack(RF, InData%RegDeterMethod) call RegPack(RF, InData%RegFunction) + call RegPack(RF, InData%RegFunctionPart) call RegPack(RF, InData%WakeRegMethod) call RegPack(RF, InData%WakeRegParam) call RegPack(RF, InData%WingRegParam) @@ -4322,6 +4329,7 @@ subroutine FVW_UnPackInputFile(RF, OutData) call RegUnpack(RF, OutData%CoreSpreadEddyVisc); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%RegDeterMethod); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%RegFunction); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%RegFunctionPart); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%WakeRegMethod); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%WakeRegParam); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%WingRegParam); if (RegCheckErr(RF, RoutineName)) return diff --git a/modules/aerodyn/src/FVW_VortexTools.f90 b/modules/aerodyn/src/FVW_VortexTools.f90 index a8589bebce..18073d8901 100644 --- a/modules/aerodyn/src/FVW_VortexTools.f90 +++ b/modules/aerodyn/src/FVW_VortexTools.f90 @@ -1448,7 +1448,7 @@ end subroutine cut_tree_segment_parallel ! --- Velocity computation ! -------------------------------------------------------------------------------- subroutine ui_tree_part(Tree, icp_end, CPs, BranchFactor, DistanceDirect, Uind, ErrStat, ErrMsg) - use FVW_BiotSavart, only: ui_part_nograd_11 + use FVW_BiotSavart, only: ui_part_nograd_11, PartRegFloorFactor type(T_Tree), target, intent(inout) :: Tree !< integer, intent(in ) :: icp_end !< Number of CPs to use