OLAF improvements for WEIS

cmake/OpenfastFortranOptions.cmake

@@ -110,7 +110,7 @@ macro(set_fast_gfortran)
   # Disable stack reuse within routines: issues seen with gfortran 9.x, but others may also exhibit
   #   see section 3.16 of https://gcc.gnu.org/onlinedocs/gcc-9.2.0/gcc.pdf
   #   and https://github.com/OpenFAST/openfast/pull/595
-  set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fstack-reuse='none'")
+  set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fstack-reuse=none")
 
   # Deal with Double/Single precision
   if (DOUBLE_PRECISION)
@@ -120,7 +120,7 @@ macro(set_fast_gfortran)
 
   # debug flags
   if(CMAKE_BUILD_TYPE MATCHES Debug)
-    set( CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} -fcheck=all -pedantic -fbacktrace " )
+    set( CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} -fcheck=all,no-array-temps -pedantic -fbacktrace " )
   endif()
 
   if(CYGWIN)

docs/conf.py

@@ -245,7 +245,10 @@ def runDoxygen(sourcfile, doxyfileIn, doxyfileOut):
 ]
 
 def setup(app):
-    app.add_css_file('css/math_eq.css')
+    try:
+        app.add_css_file('css/math_eq.css')
+    except:
+        pass
     app.add_object_type(
         "confval",
         "confval",

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

@@ -1,7 +1,7 @@
 --------------------------- OLAF (cOnvecting LAgrangian Filaments) INPUT FILE -----------------
 Free wake input file for the Helix test case
 --------------------------- GENERAL OPTIONS ---------------------------------------------------
-5       IntMethod          Integration method {5: Forward Euler 1st order, default: 5} (switch)
+5       IntMethod          Integration method {1: Runge-Kutta 4th order, 5: Forward Euler 1st order, default: 5} (switch)
 0.2     DTfvw              Time interval for wake propagation. {default: dtaero} (s)
 5       FreeWakeStart      Time when wake is free. (-) value = always free. {default: 0.0} (s)
 2.0     FullCircStart      Time at which full circulation is reached. {default: 0.0} (s)
@@ -20,11 +20,11 @@ default FreeWakeLength     Wake length that is free [integer] (number of time st
 False   FWShedVorticity    Include shed vorticity in the far wake {default: false}
 ------------------- WAKE REGULARIZATIONS AND DIFFUSION -----------------------------------------
 0       DiffusionMethod    Diffusion method to account for viscous effects {0: None, 1: Core Spreading, "default": 0}         
-0       RegDeterMethod     Method to determine the regularization parameters {0:  Manual, 1: Optimized, default: 0 }          
+0       RegDeterMethod     Method to determine the regularization parameters {0: Constant, 1: Optimized, 2: Chord-scaled, 3: dr-scaled, default: 0 }          
 2       RegFunction        Viscous diffusion function {0: None, 1: Rankine, 2: LambOseen, 3: Vatistas, 4: Denominator, "default": 3} (switch)
 0       WakeRegMethod      Wake regularization method {1: Constant, 2: Stretching, 3: Age, default: 1} (switch)
-2.0     WakeRegFactor      Wake regularization factor (m)
-2.0     WingRegFactor      Wing regularization factor (m)
+2.0     WakeRegFactor      Wake regularization factor (m or -)
+2.0     BladeRegFactor     Blade regularization factor (m or -)
 100     CoreSpreadEddyVisc Eddy viscosity in core spreading methods, typical values 1-1000 
 ------------------- WAKE TREATMENT OPTIONS ---------------------------------------------------
 False   TwrShadowOnWake    Include tower flow disturbance effects on wake convection {default:false} [only if TwrPotent or TwrShadow]
@@ -39,4 +39,7 @@ False   TwrShadowOnWake    Include tower flow disturbance effects on wake convec
 1       nVTKBlades         Number of blades for which VTK files are exported {0: No VTK per blade, n: VTK for blade 1 to n} (-) 
 2       VTKCoord           Coordinate system used for VTK export. {1: Global, 2: Hub, "default": 1} 
 1       VTK_fps            Frame rate for VTK output (frames per second) {"all" for all glue code timesteps, "default" for all OLAF timesteps} [used only if WrVTK=1] 
+0       nGridOut           Number of grid outputs
+GridName    DTOut     XStart    XEnd   nX    YStart   YEnd    nY    ZStart   ZEnd   nZ
+(-)         (s)        (m)      (m)    (-)    (m)     (m)     (-)    (m)     (m)    (-)
 ------------------------------------------------------------------------------------------------

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

@@ -30,7 +30,8 @@ General Options
 convect the Lagrangian markers. There are four options: 1) fourth-order
 Runge-Kutta *[1]*, 2) fourth-order Adams-Bashforth *[2]*, 3) fourth-order
 Adams-Bashforth-Moulton *[3]*, and 4) first-order forward Euler *[5]*. The
-default option is *[5]*. These methods are specified in :numref:`sec:vortconv`.
+default option is *[5]*. Currently only options *[1]* and *[5]* are implemented.
+These methods are specified in :numref:`sec:vortconv`.
 
 **DTfvw** [sec] specifies the time interval at which the module will update the
 wake. The time interval must be a multiple of the time step used by
@@ -111,10 +112,40 @@ for viscous diffusion. There are two options: 1) no diffusion *[0]* and 2) the
 core-spreading method *[1]*. The default option is *[0]*.
 
 **RegDetMethod** [switch] specifies which method is used to determine the
-regularization parameters. There are two options: 1) manual *[0]* and 2)
-optimized *[1]*. The manual option requires the user to specify the parameters
-listed in this subsection. The optimized option determines the parameters for
-the user.  The default option is *[0]*.
+regularization parameters. There are four options: 1) constant *[0]* and 2)
+optimized *[1]*, 3) chord *[2]*, and 4) span *[3]*. 
+The optimized option determines all the parameters in this section for the user.
+The optimized option is still work in progress and not recommended.
+The constant option requires the user to specify all the parameters present in this section.
+The default option is *[0]*.
+When **RegDetMethod==0**, the regularization parameters is set constant:
+
+.. math::
+
+   r_{c,\text{wake}}(r) = \text{WakeRegParam} 
+   ,\quad
+   r_{c,\text{blade}}(r) = \text{BladeRegParam} 
+
+When **RegDetMethod==2**, the regularization parameters is set according to the local chord:
+
+.. math::
+
+   r_{c,\text{wake}}(r) = \text{WakeRegParam} \cdot c(r)
+   ,\quad
+   r_{c,,\text{blade}}(r) = \text{BladeRegParam} \cdot c(r)
+
+When **RegDetMethod==3**, the regularization parameters is set according to the spanwise discretization:
+
+.. math::
+
+   r_{c,\text{wake}}(r) = \text{WakeRegParam} \cdot \Delta  r(r)
+   ,\quad
+   r_{c,,\text{blade}}(r) = \text{BladeRegParam} \cdot \Delta r(r)
+
+where :math:`Delta r` is the length of the spanwise station.
+
+
+
 
 **RegFunction** [switch] specifies the regularization function used to remove
 the singularity of the vortex elements, as specified in
@@ -128,11 +159,11 @@ radius (i.e., the regularization parameter). There are three options: 1)
 constant *[1]*, 2) stretching *[2]*, and 3) age *[3]*. The methods are
 described in :numref:`sec:corerad`. The default option is *[1]*.
 
-**WakeRegParam** [m] specifies the wake regularization parameter, which is the
+**WakeRegParam** [m, or -] specifies the wake regularization parameter, which is the
 regularization value used at the initialization of a vortex element. If the
 regularization method is “constant”, this value is used throughout the wake.
 
-**BladeRegParam** [m] specifies the bound vorticity regularization parameter,
+**BladeRegParam** [m, or -] specifies the bound vorticity regularization parameter,
 which is the regularization value used for the vorticity elements bound to the
 blades.
 
@@ -200,6 +231,29 @@ provided value is rounded to the nearest allowable multiple of the time step.
 The default value is :math:`1/dt_\text{fvw}`. Specifying *VTK_fps* = *[all]*,
 is equivalent to using the value :math:`1/dt_\text{aero}`.
 
+
+**nGridOut** [-] specifies the number of grid outputs. The default value is 0.
+The grid outputs are velocity fields that are exported on a regular Cartesian grid. 
+The are defined using a table that follows on the subsequent lines, with two lines of headers. 
+The user needs to specify a **GridName**, used for the VTK output filename, a time interval 
+**DTOut**, and the grid extent in each directions, e.g. **XStart**, **XEnd**, **nX**. 
+With these options, it is possible to export the velocity field at a point (**nX=nY=nZ=1**),
+a line, a plane, or a box. When the variable **DTOut** is set to "all", the AeroDyn time step is used,  when it is set to "default", the OLAF time step is used. 
+An example of input is given below: 
+
+.. code::
+
+    3       nGridOut           Number of grid outputs
+    GridName    DTOut     XStart    XEnd   nX    YStart   YEnd    nY    ZStart   ZEnd   nZ
+    (-)         (s)        (m)      (m)    (-)    (m)     (m)     (-)    (m)     (m)    (-)
+    "box"       all        -200     1000.    5    -150.   150.    20      5.     300.    30
+    "vert"      default    -200     1000.   100     0.     0.     1       5.     300.    30
+    "hori"      2.0        -200     1000.   100   -150.   150.    20     100.    100.    1
+
+In this example, the first grid, named "box", is exported at the AeroDyn time step, and consists 
+of a box of shape 5x20x30 and dimension 1200x300x295.  The two other grids are vertical and horizontal planes.
+
+
 AeroDyn15 Input File
 --------------------
 Input file modifications

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

@@ -309,9 +309,10 @@ markers are the end points of the vortex filaments. The Lagrangian convection of
 the filaments stretches the filaments and thus automatically accounts for strain
 in the vorticity equation.
 
-At present, a first-order forward Euler method is used to numerically solve the
-left-hand side of Eq. :eq:`VortFilCart` for the vortex filament location
-(**IntMethod=[5]**). This is an explicit method solved using
+At present, the Runge-Kutta 4th order (**IntMethod=[1]**) or first order forward Euler
+(**IntMethod=[5]**) methods are implemented to numerically solve the
+left-hand side of Eq. :eq:`VortFilCart` for the vortex filament location.
+In the case of the first order Euler method, the convection is then simply:
 Eq. :eq:`eq:Euler`. 
 
 .. math::
@@ -341,14 +342,6 @@ where :math:`d\psi/dt=\Omega` and
 :math:`\vec{r}(\psi,\zeta)` is the position vector of a Lagrangian
 marker, and :math:`\vec{V}[\vec{r}(\psi,\zeta)]` is the velocity.
 
-..
-   At present, first-order forward Euler method is used to numerically solve the
-   left-hand side of Eq. :eq:`VortFil_expanded` for the vortex-filament location
-   [**IntMethod=5**]. This is an explicit method solved using Eq. :eq:`Euler`.
-
-   .. math::
-      \vec{r}(\psi+\Delta\psi_i,\zeta+\Delta\zeta)  = \vec{r}(\psi,\zeta) + \vec{V}(\psi,\zeta) \Delta t
-      :label: Euler
 
 Induced Velocity and Velocity Field
 -----------------------------------
@@ -577,6 +570,8 @@ Here, :math:`\epsilon` is the vortex-filament strain, :math:`l` is the filament
 length, and :math:`\Delta l` is the change of length between two time steps. The
 integral in Eq. :eq:`stretch` represents strain effects.
 
+This option is not yet implemented.
+
 Wake Age / Core-Spreading
 ^^^^^^^^^^^^^^^^^^^^^^^^^
 
@@ -599,6 +594,16 @@ vorticity itself or between the wake vorticity and the background flow. It is
 often referred to as the core-spreading method. Setting **DiffusionMethod=[1]**
 is the same as using the wake age method (**WakeRegMethod=[3]**).
 
+The time evolution of the core radius is implemented as:
+
+.. math::
+
+    \frac{d r_c}{dt} = \frac{2\alpha\delta\nu}{r_c(t)}
+
+and :math:`\frac{d r_c}{dt}=0` on the blades.
+
+
+
 Stretching and Wake Age
 ^^^^^^^^^^^^^^^^^^^^^^^
 
@@ -610,6 +615,8 @@ Eq. :eq:`stretchandage`.
    r_c(\zeta,\epsilon) = \sqrt{r_{c0}^2 + 4\alpha\delta\nu \zeta  \big(1+\epsilon\big)^{-1} }
    :label: stretchandage
 
+This option is not yet implemented.
+
 .. _sec:diffusion:
 
 Diffusion

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

@@ -11,6 +11,9 @@ parameter is available in the OLAF input file, in which case the VTK files are
 written to the folder ``vtk_fvw``, or the primary ``.fst`` file, in which case
 the VTK files are written to the folder ``vtk``.
 
+Velocity field outputs can be exported as VTK files. The user can control these
+outputs using **nGridOut** and the subsequent table.
+
 
 Results File
 ------------

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

@@ -1,7 +1,12 @@
+
+Working with OLAF
+=================
+
+
 .. _Running-OLAF:
 
 Running OLAF
-============
+~~~~~~~~~~~~
 
 As OLAF is a module of OpenFAST, the process of downloading, compiling,
 and running OLAF is the same as that for OpenFAST. Such instructions are
@@ -10,3 +15,83 @@ available in the :ref:`installation` documentation.
 .. note::
    To improve the speed of FVW module, the user may wish to compile with
    `OpenMP`.  To do so, add the `-DOPENMP=ON` option with CMake.
+
+
+Guidelines
+~~~~~~~~~~
+
+Most options of OLAF can be left to default. The results will depend on the time discretization, wake length, and regularization parameters. We provide guidelines for these parameters in this section, together with a simple python code to compute these parameters.
+Please check this section again as we might further refine our guidelines with time.
+
+
+**Time step and wake length**
+We recommend to set OLAF's time step (**DTfvw**) such that it corresponds to :math:`\Delta \psi = 6` degrees of one rotor revolution:
+
+.. math::
+
+    \Delta t
+    = \frac{\Delta \psi_\text{rad}}{\Omega_\text{rad/s}}
+    = \frac{\Delta \psi_\text{deg}}{6 \times \Omega_\text{RPM}}
+
+If the structural solver requires a smaller time step, the time step for the glue code can be set to a different value than **DTfvw** as long as **DTfvw** is a multiple of the glue code time step.
+
+
+We recommend to set the near wake length to the number of time steps needed to reach two rotor revolutions. For the far wake, we recommend 10 rotor revolutions. 
+For the free far-wake, we recommend to set the distance to a value somewhere between 25% and 50% of the far wake length, (e.g. 3 revolutions).
+
+The following python script computes the parameters according to these guidelines.
+
+.. code::
+
+   def OLAFParams(omega_rpm, deltaPsiDeg=6, nNWrot=2, nFWrot=10, nFWrotFree=3, nPerAzimuth=None):
+       """ 
+       Computes recommended time step and wake length based on the rotational speed in RPM
+
+       INPUTS:
+        - omega_rpm: rotational speed in RPM
+        - deltaPsiDeg : azimuthal discretization in deg
+        - nNWrot : number of near wake rotations
+        - nFWrot : total number of far wake rotations
+        - nFWrotFree : number of far wake rotations that are free
+
+           deltaPsiDeg  -  nPerAzimuth
+                5            72    
+                6            60    
+                7            51.5  
+                8            45    
+       """
+       omega = omega_rpm*2*np.pi/60
+       T = 2*np.pi/omega
+       if nPerAzimuth is not None:
+           dt_wanted    = np.around(T/nPerAzimuth,3)
+       else:
+           dt_wanted    = np.around(deltaPsiDeg/(6*omega_rpm),3)
+           nPerAzimuth = int(2*np.pi /(deltaPsiDeg*np.pi/180))
+
+       nNWPanel     = nNWrot*nPerAzimuth
+       nFWPanel     = nFWrot*nPerAzimuth
+       nFWPanelFree = nFWrotFree*nPerAzimuth
+
+       print(dt_wanted              , '  DTfvw')
+       print(int      (nNWPanel    ), '  nNWPanel      ')
+       print(int      (nFWPanel    ), '  WakeLength    ')
+       print(int      (nFWPanelFree), '  FreeWakeLength')
+
+       return dt_wanted, nNWPanel, nFWPanel, nFWPanelFree
+
+
+**Regularization parameters**
+
+One critical parameter of vortex methods is the regularization parameter, also referred to as core radius. We currently recommend to set the regularization parameter as a fraction of the spanwise discretization, that is: **RegDetMethod=3** , **WakeRegFactor=0.6**, **WingRegFactor=0.6**.
+We will likely update these guidelines in the future.
+
+
+We also recommend to have the regularization increasing with downstream distance:
+**WakeRegMethod=3**. 
+
+The factor with which the regularization parameter will increase with downstream distance can be set as
+**CoreSpreadEddyVisc=1000** for modern multi-MW turbines. Further guidelines will follow for this parameter in the future. 
+
+
+
+