[OLAF] Implementation of a source panel method

docs/source/acknowledgements.rst

@@ -11,8 +11,37 @@ from U.S. Department of Energy Wind Energy Technology Office through the
 NREL gratefully acknowledges development contributions from the following
 organizations:
 
-- Envision Energy USA, Ltd
+- Accelerate Wind
+- AirLoom Energy
+- ARPA-E
 - Brigham Young University
+- CENER
+- Delft University of Technology
+- DOE Water Power Technologies Office (WPTO)
+- Equinor
+- Envision Energy USA, Ltd
+- Institute for Energy Technology
+- GoogleX
+- Hayman Consulting
+- Norwegian University of Science and Technology
+- NOWRDC
+- Octue
+- Orcina
+- Oregon State University
+- Principle Power
+- RRD Engineering
+- Sandia National Laboratory
+- SENER
+- Shell
+- SOWENTO
+- Technical University of Denmark
+- TODA Corporation
+- TotalEnergies
+- University of Strathclyde
+- University of Stuttgart
+- University of Texas-Austin
+- Vestas
+- Windward Engineering
 
 NREL gratefully acknowledges additional development support through designation
 as an `Intel® Parallel Computing Center (IPCC) <https://software.intel.com/en-us/ipcc>`_.

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

@@ -43,4 +43,8 @@ default VTK_fps            - Frame rate for VTK output (frames per second) {"all
 0       nGridOut           - Number of grid outputs
 GridName  GridType  TStart  TEnd     DTGrid    XStart    XEnd   nX    YStart   YEnd    nY    ZStart   ZEnd   nZ
 (-)         (-)      (s)     (s)      (s)        (m)      (m)    (-)    (m)     (m)     (-)    (m)     (m)    (-)
-------------------------------------------------------------------------------------------------
+===============================================================================================
+--------------------------- ADVANCED OPTIONS --------------------------------------------------
+===============================================================================================
+! Advanced options may be placed here in arbitrary order, with the regular format:
+! Value   Key             - Comment

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

@@ -312,6 +312,118 @@ of a box of shape 5x20x30 and dimension 1200x300x295.  The grid contains both th
 The two other grids are vertical and horizontal planes containing only the velocity.
 
 
+Advanced Options
+~~~~~~~~~~~~~~~~
+
+
+Advanced options (typically used for developers or beta features) can be placed at the end of the OLAF input file:
+
+
+- These options use the regular format:  `Value Key - Comment`.
+- They can be placed in arbitrary order.
+- They can be commented out with a `!` character at the beginning of the line.
+- Blank lines or unsupported options are ignored (display to screen).
+- If not provided, a default value will be used, but the "DEFAULT" keyword is not supported for the advanced options.
+- These are **beta features** and should mostly be used by developers.
+
+
+The end of the OLAF input file would look as follows:
+
+.. code::
+
+   [...]
+   GridName  GridType  TStart  TEnd     DTGrid    XStart    XEnd   nX    YStart   YEnd    nY    ZStart   ZEnd   nZ
+   (-)         (-)      (s)     (s)      (s)        (m)      (m)    (-)    (m)     (m)     (-)    (m)     (m)    (-)
+   ===============================================================================================
+   --------------------------- ADVANCED OPTIONS --------------------------------------------------
+   ===============================================================================================
+   ! Advanced options may be placed here in arbitrary order, with the regular format:
+   ! Value1   Key1            - Comment1
+   ! Value2   Key2            - Comment2
+   ! Lines starting with `!` are ignored, empty lines are ignored
+   [...] etc
+
+Currently, the advanced options supported are as follows:
+
+.. code::
+
+   ===============================================================================================
+   --------------------------- ADVANCED OPTIONS --------------------------------------------------
+   ===============================================================================================
+   "Panels.vtk"  SrcPnlFile     - Name of VTK file containing source panels {default: ""}
+   1             nSrcPnlUpdate  - How often do src panel updates (in time steps of OLAF), {default: 1}
+   True          Induction      - Compute induced velocities from wake to blade and wake to wake, {default: True} 
+   True          InductionAtCP  - Compute induced velocities at nodes or CP, {default: True} 
+   True          WakeAtTE       - Start the wake at the trailing edge, or at the LL, {default: True} 
+   False         DStallOnWake   - Dynamic stall has influence on wake, {default: False} 
+   0.75          kFrozenNWStart - Fraction of wake induced velocity at start of frozen wake, {default: 0.75} 
+   0.5           kFrozenNWEnd   - Fraction of wake induced velocity at end of frozen wake, {default: 0.5} 
+   0.0           zGround        - Ground height, used to enforce that no vortices go into the ground {default: 0.0}
+   0.1           zGroundPush    - Ground push back, vortices that are lower than zGround are placed back at zGroundPush {default: 0.1}
+
+
+
+**SrcPanlFile** [string] specifies the name of the VTK file used for the source panel method.
+The VTK file should be a legacy ASCII VTK file, with DATASET POLYDATA POINTS and POLYGONS.
+A sample VTK file is provided below for two panels forming a regular grid in the XY plane. 
+The connectivity of the polygons needs to be such that the normal points away from the body and
+into the fluid. 
+In the example below, the polygons are defined clockwise when viewed from above, which results in the normals internally defined by OLAF pointing in the `+z` direction. This configuration is typical of a bottom wall with fluid above it: the OLAF normals point from the wall toward the fluid. However, when applying the right-hand-rule convention, the resulting normals point from the fluid toward the interior of the body.
+When using **WrVTK**, OLAF will write separate VTK files containing various information related to the source panels, such as the pressure coefficient, area, force per area, normals. Looking at the orientation of the normals is extremely important (In Paraview, select 3D Glyph, Arrows, and select the `Normals` output, scaled by the `Normals`). 
+The BodyID CELL_DATA can be used to separate different patches, which can help the post processing. This is optional, OLAF doesn't use it currently, but it is written back in the output files of OLAF.
+
+Curious readers may look at the unittest `Test_SrcPnl_Sphere` of OLAF that test the pressure coefficient
+about a sphere.
+
+.. code::
+
+   # vtk DataFile Version 2.0
+   Comment
+   ASCII
+   DATASET POLYDATA
+   POINTS 6 double
+   0.0 0.0 0.0
+   0.0 25.0 0.0
+   0.0 50.0 0.0
+   10.0 0.0 0.0
+   10.0 25.0 0.0
+   10.0 50.0 0.0
+
+   POLYGONS 2 10
+   4 0 1 4 3
+   4 1 2 5 4
+
+   CELL_DATA 2
+   SCALARS BodyID int
+   LOOKUP_TABLE default
+   0
+   0
+
+
+
+**nSrcPnlUpdate** [int] Define how often the source panel updates (in time steps of OLAF), the default is `1`, at each time step.
+
+**Induction**  [switch] Compute induced velocities, otherwise all induced velocity are 0 (no wake, no panels, etc) ! Default is `True`.
+
+**InductionAtCP** [switch] When performing the lifting line calculations, compute induced velocities at nodes (False) or at control points (CP, True). Default is `True`.
+
+**WakeAtTE** [switch] Start the wake at the trailing edge (True), or directly at the LL (False, no chordwise panel). Default is `True`.
+
+**DStallOnWake** [switch] Include the influence of the dynamic stall on the wake (True), i.e. use the dynamic Cl to update the circulation of the lifting line and wake panel. Default is `False`.
+
+**kFrozenNWStart** [float] Fraction of wake induced velocity at start of frozen wake, default is `0.75`. See OLAF theory related to Frozen NW, :numref:`sec:vortconvfrozen`.
+
+**kFrozenNWEnd** [float] Fraction of wake induced velocity at end of frozen wake, default is `0.5`. See OLAF theory related to Frozen NW, :numref:`sec:vortconvfrozen`.
+
+**zGround** [float] Height in meters below which vortex points are not allowed to be and if any are present they will be pushed back above the ground at the height defined by **zGroundPush** (in meters). Default is `0.0`.
+For MHK, the sea bed location, based on the water depth is added to **zGround**.
+
+**zGroundPush** [float] Ground push back, see **zGround**. Default is `0.1`.
+
+
+
+
+
 AeroDyn Input File
 --------------------
 Input file modifications

modules/aerodyn/src/AeroDyn.f90

@@ -412,7 +412,12 @@ subroutine AD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut
 
 
       ! set the rest of the parameters
-   p%Skew_Mod = InputFileData%Skew_Mod
+   ! NOTE: some parameters need to be set even if no rotors
+   p%Skew_Mod      = InputFileData%Skew_Mod
+   p%UA_Flag       = InputFileData%UA_Init%UAMod > UA_None
+   p%CompAeroMaps  = InitInp%CompAeroMaps
+   p%DT            = InputFileData%DTAero
+   p%Wake_Mod      = InputFileData%Wake_Mod
    do iR = 1, nRotors
       p%rotors(iR)%AeroProjMod = AeroProjMod(iR)
       call WrScr('   AeroDyn: projMod: '//trim(num2lstr(p%rotors(iR)%AeroProjMod)))
@@ -1472,9 +1477,8 @@ subroutine SetParameters( InitInp, InputFileData, RotData, p, p_AD, ErrStat, Err
 
    p%MHK              = InitInp%MHK
 
-   p_AD%DT            = InputFileData%DTAero
-   p_AD%Wake_Mod      = InputFileData%Wake_Mod
    p%DBEMT_Mod        = InputFileData%DBEMT_Mod
+
    p%TwrPotent        = InputFileData%TwrPotent
    p%TwrShadow        = InputFileData%TwrShadow
    p%TwrAero          = InputFileData%TwrAero
@@ -1765,9 +1769,11 @@ subroutine AD_End( u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg )
       if (p%Wake_Mod == WakeMod_FVW ) then
 
          if ( p%UA_Flag ) then
-            do iW=1,p%FVW%nWings
-               call UA_End(m%FVW%W(iW)%p_UA)
-            enddo
+            if (allocated(m%FVW%W)) then
+               do iW=1,p%FVW%nWings
+                  call UA_End(m%FVW%W(iW)%p_UA)
+               enddo
+            endif
          end if
 
          call FVW_End( m%FVW_u, p%FVW, x%FVW, xd%FVW, z%FVW, OtherState%FVW, m%FVW_y, m%FVW, ErrStat, ErrMsg )
@@ -5158,11 +5164,25 @@ SUBROUTINE Init_OLAF( InputFileData, u_AD, u, p, x, xd, z, OtherState, m, ErrSta
    allocate(InitInp%W(nWings)        , STAT = ErrStat2); ErrMsg2='Allocate W'; if(Failed()) return
    allocate(InitInp%WingsMesh(nWings), STAT = ErrStat2); ErrMsg2='Allocate Wings Mesh'; if(Failed()) return
 
+   ! --- Default inputs (not per retor) 
+   ! UA and inputs that are not per rotor
+   InitInp%UA_Flag    = p%UA_Flag
+   ! Important inputs if no rotors are present!
+   if (size(p%rotors)==0) then
+      InitInp%numBladeNodes = 0
+      InitInp%AirDens       = InputFileData%AirDens
+      InitInp%KinVisc       = InputFileData%KinVisc
+      InitInp%MHK           = 0                                    ! TODO this is an initinp of AeroDyn
+      InitInp%WtrDpth       = 0.0_ReKi                             ! TODO this is an initinp of AeroDyn
+      InitInp%RootName      = p%RootName(1:len_trim(p%RootName)-2) ! Removing "AD"
+   endif
+
    ! --- Inputs per wings/blades
    iW_incr=0
    do iR=1, size(p%rotors)
 
       InitInp%numBladeNodes  = p%rotors(iR)%numBlNds ! TODO TODO TODO per wing
+      InitInp%AirDens        = p%rotors(iR)%AirDens
       InitInp%KinVisc        = p%rotors(iR)%KinVisc
       InitInp%MHK            = p%rotors(iR)%MHK
       InitInp%WtrDpth        = p%rotors(iR)%WtrDpth
@@ -5225,9 +5245,6 @@ SUBROUTINE Init_OLAF( InputFileData, u_AD, u, p, x, xd, z, OtherState, m, ErrSta
          if(Failed()) return
 
       enddo ! iB, blades
-
-      ! Unsteady Aero Data
-      InitInp%UA_Flag    = p%UA_Flag
       call UA_CopyInitInput(InputFileData%UA_Init, InitInp%UA_Init, MESH_NEWCOPY, ErrStat2, ErrMsg2)
 
       iW_incr = iW_incr+p%rotors(iR)%numBlades

modules/aerodyn/src/AeroDyn_Driver_Subs.f90

@@ -541,7 +541,7 @@ subroutine Init_ADI_ForDriver(iCase, ADI, dvr, FED, dt, needInitIW, errStat, err
       InitInp%AD%defPvap     = dvr%Pvap
       InitInp%AD%WtrDpth     = dvr%WtrDpth
       InitInp%AD%MSL2SWL     = dvr%MSL2SWL
-      InitInp%AD%CompSeaSt   = dvr%SS_InitInp%CompSeaSt
+      InitInp%AD%CompSeaSt   = dvr%SS_InitInp%CompSeaSt /=0
       ! Init data per rotor
       allocate(InitInp%AD%rotors(dvr%numTurbines), stat=errStat) 
       if (errStat/=0) then
@@ -1038,7 +1038,7 @@ subroutine Dvr_ReadInputFile(fileName, dvr, errStat, errMsg )
    call ParseCom(FileInfo_In, CurLine, Line, errStat2, errMsg2, unEc); if (Failed()) return
    call ParseVar(FileInfo_In, CurLine, "compInflow", dvr%IW_InitInp%compInflow  , errStat2, errMsg2, unEc); if (Failed()) return
    call ParseVar(FileInfo_In, CurLine, "InflowFile", dvr%IW_InitInp%InputFile, errStat2, errMsg2, unEc, IsPath=.true.); if (Failed()) return
-   if (dvr%IW_InitInp%compInflow==0) then
+   if (dvr%IW_InitInp%compInflow/=1) then
       call ParseVar(FileInfo_In, CurLine, "HWindSpeed", dvr%IW_InitInp%HWindSpeed  , errStat2, errMsg2, unEc); if (Failed()) return
       call ParseVar(FileInfo_In, CurLine, "RefHt"     , dvr%IW_InitInp%RefHt       , errStat2, errMsg2, unEc); if (Failed()) return
       call ParseVar(FileInfo_In, CurLine, "PLExp"     , dvr%IW_InitInp%PLExp       , errStat2, errMsg2, unEc); if (Failed()) return
@@ -1544,6 +1544,9 @@ subroutine Dvr_InitializeOutputs(nWT, out, numSteps, errStat, errMsg)
             end do ! i
             write (out%unOutFile(iWT),'()')
          enddo
+         if (nWT==0) then
+            ! Special case, no turbine, TODO
+         endif
       endif
 
       ! --- Binary
@@ -1745,7 +1748,11 @@ subroutine Dvr_WriteOutputs(nt, t, dvr, out, yADI, SeaSt, errStat, errMsg)
    END IF
 
    ! Packing all outputs excpet time into one array
-   nAD = size(yADI%AD%rotors(1)%WriteOutput)
+   if (size(yADI%AD%rotors)>=1) then
+      nAD = size(yADI%AD%rotors(1)%WriteOutput)
+   else
+      nAD = 0
+   endif
    nIW = size(yADI%IW_WriteOutput)
    nSS = size(SeaSt%y%WriteOutput)
    nDV = out%nDvrOutputs
@@ -1873,11 +1880,15 @@ subroutine setVTKParameters(DVR_Outs, dvr, ADI, errStat, errMsg, dirname)
    ! Get radius for ground (blade length + hub radius):
    GroundRad = MaxBladeLength + MaxTwrLength+ DVR_Outs%VTKHubRad
    ! write the ground or seabed reference polygon:
-   RefPoint(1:2) = dvr%WT(1)%originInit(1:2)
-   do iWT=2,dvr%numTurbines
-      RefPoint(1:2) = RefPoint(1:2) + dvr%WT(iWT)%originInit(1:2)
-   end do
-   RefPoint(1:2) = RefPoint(1:2) / dvr%numTurbines
+   if (dvr%numTurbines>0) then
+      RefPoint(1:2) = dvr%WT(1)%originInit(1:2)
+      do iWT=2,dvr%numTurbines
+         RefPoint(1:2) = RefPoint(1:2) + dvr%WT(iWT)%originInit(1:2)
+      end do
+      RefPoint(1:2) = RefPoint(1:2) / dvr%numTurbines
+   else
+      RefPoint(1:3) = 0.0_ReKi
+   endif
 
    RefPoint(3) = 0.0_ReKi
    RefLengths  = GroundRad  + sqrt((WorldBoxMax(1)-WorldBoxMin(1))**2 + (WorldBoxMax(2)-WorldBoxMin(2))**2)

modules/aerodyn/src/AeroDyn_Inflow.f90

@@ -86,12 +86,18 @@ subroutine ADI_Init(InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut
    InitOut%Ver = InitOut_AD%ver
    ! Add writeoutput units and headers to driver, same for all cases and rotors!
    !TODO: this header is too short if we add more rotors.  Should also add a rotor identifier
-   call concatOutputHeaders(InitOut%WriteOutputHdr, InitOut%WriteOutputUnt, InitOut_AD%rotors(1)%WriteOutputHdr, InitOut_AD%rotors(1)%WriteOutputUnt, errStat2, errMsg2); if(Failed()) return
+   if (size(InitOut_AD%rotors)>=1) then
+      call concatOutputHeaders(InitOut%WriteOutputHdr, InitOut%WriteOutputUnt, InitOut_AD%rotors(1)%WriteOutputHdr, InitOut_AD%rotors(1)%WriteOutputUnt, errStat2, errMsg2); if(Failed()) return
+   endif
 
    ! --- Initialize grouped outputs
-   !TODO: assumes one rotor
-   p%NumOuts = p%AD%rotors(1)%NumOuts + p%AD%rotors(1)%BldNd_TotNumOuts + m%IW%p%NumOuts
-   call AllocAry(y%WriteOutput, p%NumOuts, 'WriteOutput', errStat2, errMsg2); if (Failed()) return
+   if (size(InitOut_AD%rotors)>=1) then
+      !TODO: assumes one rotor
+      p%NumOuts = p%AD%rotors(1)%NumOuts + p%AD%rotors(1)%BldNd_TotNumOuts + m%IW%p%NumOuts
+      call AllocAry(y%WriteOutput, p%NumOuts, 'WriteOutput', errStat2, errMsg2); if (Failed()) return
+   else
+      p%NumOuts = m%IW%p%NumOuts
+   endif
 
    ! --- Initialize outputs
    call AllocAry(y%IW_WriteOutput, size(m%IW%y%WriteOutput),'IW_WriteOutput', errStat2, errMsg2); if(Failed()) return
@@ -297,21 +303,26 @@ subroutine ADI_CalcOutput(t, u, p, x, xd, z, OtherState, y, m, errStat, errMsg)
    y%PLExp = m%IW%PLExp
 
    ! --- Set outputs
-   !TODO: this assumes one rotor!!!
-   associate(AD_NumOuts => p%AD%rotors(1)%NumOuts + p%AD%rotors(1)%BldNd_TotNumOuts, &
-             IW_NumOuts => m%IW%p%NumOuts)
-      y%WriteOutput(1:IW_NumOuts) = y%IW_WriteOutput(1:IW_NumOuts)
-      y%WriteOutput(IW_NumOuts+1:p%NumOuts) = y%AD%rotors(1)%WriteOutput(1:AD_NumOuts)
-   end associate
-
-   !----------------------------------------------------------------------------
-   ! Store hub height velocity calculated in CalcOutput
-   !----------------------------------------------------------------------------
+   if (size(p%AD%rotors)>=1) then
+      !TODO: this assumes one rotor!!!
+      associate(AD_NumOuts => p%AD%rotors(1)%NumOuts + p%AD%rotors(1)%BldNd_TotNumOuts, &
+                IW_NumOuts => m%IW%p%NumOuts)
+         y%WriteOutput(1:IW_NumOuts) = y%IW_WriteOutput(1:IW_NumOuts)
+         y%WriteOutput(IW_NumOuts+1:p%NumOuts) = y%AD%rotors(1)%WriteOutput(1:AD_NumOuts)
+      end associate
+
+      !----------------------------------------------------------------------------
+      ! Store hub height velocity calculated in CalcOutput
+      !----------------------------------------------------------------------------
+
+      if (p%storeHHVel) then
+         do iWT = 1, size(u%AD%rotors)
+            y%HHVel(:,iWT) = m%AD%Inflow(1)%RotInflow(iWT)%InflowOnHub(:,1)
+         end do
+      endif
 
-   if (p%storeHHVel) then
-      do iWT = 1, size(u%AD%rotors)
-         y%HHVel(:,iWT) = m%AD%Inflow(1)%RotInflow(iWT)%InflowOnHub(:,1)
-      end do
+   else
+      y%WriteOutput(1:p%NumOuts) = y%IW_WriteOutput(1:m%IW%p%NumOuts)
    endif
 
 contains