OpenFAST · rafmudaf · Sep 26, 2022 · Jan 29, 2022 · Jan 29, 2022 · Jan 29, 2022
diff --git a/docs/source/user/aerodyn/bibliography.bib b/docs/source/user/aerodyn/bibliography.bib
@@ -26,6 +26,18 @@ @book{ad-Branlard:book
 }
 
 
+@article{ad-Branlard:2022,
+	author = {E Branlard and B Jonkman and G R Pirrung and K Dixon and J Jonkman},
+	title = {Dynamic inflow and unsteady aerodynamics models for modal and stability analyses in {OpenFAST}},
+	doi = {10.1088/1742-6596/2265/3/032044},
+	year = 2022,
+	publisher = {{IOP} Publishing},
+	volume = {2265},
+	number = {3},
+	pages = {032044},
+	journal = {Journal of Physics: Conference Series}
+}
+
 @article{ad-Hansen:book,
     author = {Hansen, M. O. L. and S{\o}rensen, J. N. and Voutsinas, S. and S{\o}rensen, N. and Madsen, H. Aa.},
     doi = {10.1016/j.paerosci.2006.10.002},

diff --git a/docs/source/user/aerodyn/input.rst b/docs/source/user/aerodyn/input.rst
@@ -67,7 +67,7 @@ program).
 Set ``WakeMod`` to 0 if you want to disable rotor wake/induction effects or 1 to
 include these effects using the (quasi-steady) BEM theory model. When
 ``WakeMod`` is set to 2, a dynamic BEM theory model (DBEMT) is used (also
-referred to as dynamic inflow or dynamic wake model).  When ``WakeMod`` is set
+referred to as dynamic inflow or dynamic wake model, see :numref:`AD_DBEMT`).  When ``WakeMod`` is set
 to 3, the free vortex wake model is used, also referred to as OLAF (see
 :numref:`OLAF`). ``WakeMod`` cannot be set to 2 or 3 during linearization
 analyses.
@@ -179,13 +179,20 @@ Dynamic Blade-Element/Momentum Theory Options
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 The input parameters in this section are used only when ``WakeMod = 2``.
+The theory is described in :numref:`AD_DBEMT`.
+
+There are three options available for ``DBEMT_Mod``:
+
+- ``1``: discrete-time Oye's model, with constant :math:`\tau_1`
+- ``2``: discrete-time Oye's model, with varying :math:`\tau_1`, automatically adjusted based on inflow. (recommended for time-domain simulations)
+- ``3``: continuous-time Oye's model, with constant :math:`\tau_1` (recommended for linearization)
+
+For ``DBEMT_Mod=1`` or ``DBEMT_Mod=3`` it is the user responsability to set the value of :math:`\tau_1` (i.e. ``tau1_const``) according to the expression given in :numref:`AD_DBEMT`, using an estimate of what the mean axial induction (:math:`\overline{a}`) and the mean relative wind velocity across the rotor (:math:`\overline{U_0}`) are for a given simulation. 
+
+The option ``DBEMT_Mod=3`` is the only one that can be used for linearization.
+
 
-Set ``DBEMT_Mod`` to 1 for the constant-tau1 model, set ``DBEMT_Mod`` to 2
-to use a model where tau1 varies with time, or set ``DBEMT_Mod`` to 3
-to use a continuous-state model with constant tau1.
 
-If ``DBEMT_Mod=1`` (constant-tau1 model) or ``DBEMT_Mod=3`` (continuous-state constant-tau1 model), 
-set ``tau1_const`` to the time constant to use for DBEMT.
 
 OLAF -- cOnvecting LAgrangian Filaments (Free Vortex Wake) Theory Options
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -207,15 +214,15 @@ The input parameters in this section are used only when ``AFAeroMod
 
 ``UAMod`` determines the UA model. It has the following options:
 
-- ``1``: the original theoretical developments of B-L (**not currently functional**), 
-- ``2``: the extensions to B-L developed by González 
-- ``3``: the extensions to B-L developed by Minnema/Pierce
-- ``4``: a continuous-state model developed by Hansen, Gaunna, and Madsen (HGM)
-- ``5``: a model similar to HGM with an additional state for vortex generation
-- ``6``: Oye's dynamic stall model
-- ``7``: Boeing-Vertol model
+- ``1``: the discrete-time model of Beddoes-Leishman (B-L) (**not currently functional**), 
+- ``2``: the extensions to B-L developed by González  (changes in Cn, Cc, Cm)
+- ``3``: the extensions to B-L developed by Minnema/Pierce (changes in Cc and Cm)
+- ``4``: 4-states continuous-time B-L model developed by Hansen, Gaunna, and Madsen (HGM). NOTE: might require smaller time steps until a stiff integrator is implemented.
+- ``5``: 5-states continuous-time B-L model similar to HGM with an additional state for vortex generation
+- ``6``: 1-state continuous-time developed by Oye
+- ``7``: discrete-time Boeing-Vertol (BV) model
 
-The models are described in :numref:`AD_UA`.
+Linearization is supported with ``UAMod=4,5,6`` (which use continuous-time states) but not with the other models. The different models are described in :numref:`AD_UA`.
 
 
 **While all of the UA models are documented in this

diff --git a/docs/source/user/aerodyn/theory.rst b/docs/source/user/aerodyn/theory.rst
@@ -14,6 +14,98 @@ Steady BEM
 The steady blade element momentum (BEM) equations are solved as a constrained equation, and the formulation follows the description from Ning :cite:`ad-Ning:2014`.
 
 
+
+.. _AD_DBEMT:
+
+Dynamic BEM Theory (DBEMT)
+~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+
+
+Two equivalent versions of Oye's dynamic inflow model are implemented in AeroDyn.
+The first one uses discrete time, it can be used with the constant-tau1 model 
+(``DBEMT_Mod=1``) or the varying-tau1 model (``DBEMT_Mod=2``), but it cannot be used for linearization.
+The second version uses a continuous-time state-space formulation  (``DBEMT_Mod=1``), it assumes a constant-tau1, and can be used for linearization.  
+For a same value of :math:`\tau_1`, the discrete-time and continuous-time formulations returns exactly the same results.
+
+
+
+
+
+Oye's dynamic inflow model consists of two first-order differential equations (see :cite:`ad-Branlard:book`):
+
+.. math::
+   \begin{align}
+       \boldsymbol{W}_\text{int}+\tau_1    \boldsymbol{\dot{W}}_\text{int}
+           &=
+       \boldsymbol{W}_\text{qs} + k \tau_1 \boldsymbol{\dot{W}}_\text{qs} \\
+       \boldsymbol{W}+\tau_2 \boldsymbol{\dot{W}}
+           &=
+       \boldsymbol{W}_\text{int}
+   \end{align}
+
+where 
+:math:`\boldsymbol{W}` is the dynamic induction vector at the rotor (at a given blade position and radial position), 
+:math:`\boldsymbol{W}_\text{qs}` is the quasi-steady induction, 
+:math:`\boldsymbol{W}_\text{int}` is an intermediate value coupling the quasi-steady and the actual inductions (may be discontinuous if the quasi-steady indution is discontinuous).
+and
+:math:`(\dot{\ })` represents the time derivative.
+The coupling constant :math:`k`, with values between 0 and 1, is usually chosen as :math:`k=0.6`.
+Oye's dynamic inflow model relies on two time constants, :math:`\tau_1` and :math:`\tau_2` :
+
+.. math::
+        \tau_1=\frac{1.1}{1-1.3 \min(\overline{a},0.5)} \frac{R}{\overline{U}_0}
+        , \qquad
+        \tau_2 =\left[ 0.39-0.26\left(\frac{r}{R}\right)^2\right] \tau_1
+
+where :math:`R` is the rotor radius, :math:`\overline{U}_0` is the average wind speed over the rotor, :math:`\overline{a}` is the average axial induction over the rotor, and :math:`r` is the radial position along the blade.
+For ``DBEMT_Mod=1`` or ``DBEMT_Mod=3``, the user needs to provide the value of :math:`\tau_1`.
+
+
+
+
+The continuous-time state-space formulation of the dynamic inflow model (``DBEMT_Mod=3``) was derived in :cite:`ad-Branlard:2022`.
+
+.. math::
+   \begin{align}
+      \begin{bmatrix}
+      \boldsymbol{\dot{W}}_\text{red}\\
+      \boldsymbol{\dot{W}}\\
+      \end{bmatrix}
+      =
+      \begin{bmatrix}
+      -\frac{1}{\tau_1}\boldsymbol{I}_2 & \boldsymbol{0} \\
+       \frac{1}{\tau_2}\boldsymbol{I}_2 &
+      -\frac{1}{\tau_2}\boldsymbol{I}_2 \\
+      \end{bmatrix}
+      \begin{bmatrix}
+      \boldsymbol{W}_\text{red}\\
+      \boldsymbol{W}\\
+      \end{bmatrix}
+      +
+      \begin{bmatrix}
+       \frac{1-k}{\tau_1} \\
+       \frac{k}{\tau_2}\\
+      \end{bmatrix}
+     \boldsymbol{W}_\text{qs}
+   \end{align}
+
+where 
+:math:`\boldsymbol{I}_2` is the 2x2 identity matrix,
+:math:`\boldsymbol{W}_\text{red}` is the reduced induction which is a continuous, scaled, and lagged version of the quasi-steady induction, defined as:
+
+.. math::
+    \boldsymbol{W}_\text{int} = \boldsymbol{W}_\text{red} + k \boldsymbol{W}_\text{qs} 
+
+
+The discrete-time version of the model is documented in the unpublished manual of DBEMT.
+The current discrete-time formulation is complex and in the future it can be simplified by using :math:`\boldsymbol{W}_\text{red}`.
+
+
+
+
+
+
 .. _AD_twr_shadow:
 
 Tower shadow models

diff --git a/docs/source/user/aerodyn/theory_ua.rst b/docs/source/user/aerodyn/theory_ua.rst
@@ -197,8 +197,11 @@ Two variants are implemented in the Unsteady Aerodynamic module. These two (comp
 Beddoes-Leishman 4-states model (UAMod=4)
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-The 4-states (incompressible) dynamic stall model from Hansen-Gaunaa-Madsen (HGM) is described in :cite:`ad-Hansen:2004` and enabled using ``UAMod=4``.  The model uses :math:`C_l` as main physical quantity. 
-Linearization of the model will be available in the future.
+The 4-states (incompressible) dynamic stall model as implemented in OpenFAST is described in :cite:`ad-Branlard:2022` (the model differs slithgly from the original formulation from Hansen-Gaunaa-Madsen (HGM) :cite:`ad-Hansen:2004`).
+The model is enabled using ``UAMod=4``.  The model uses :math:`C_l` as main physical quantity. 
+Linearization of the model is available.
+
+NOTE: this model might require smaller time steps until a stiff integrator is implemented in AeroDyn-UA.
 
 
 **State equation:**
@@ -225,6 +228,9 @@ with
     \end{aligned}
 
 
+
+
+
 **Output equation:**
 The unsteady airfoil coefficients
 :math:`C_{l,\text{dyn}}`, :math:`C_{d,\text{dyn}}`,
@@ -233,8 +239,9 @@ The unsteady airfoil coefficients
 .. math::
 
    \begin{aligned}
-       C_{l,\text{dyn}}(t) &= x_4 (\alpha_E-\alpha_0) C_{l,\alpha} +  (1-x_4) C_{l,{fs}}(\alpha_E)+ \pi T_u \omega   \\
-       C_{d,\text{dyn}}(t) &=  C_d(\alpha_E) + (\alpha_{ac}-\alpha_E) C_{l,\text{dyn}} + \left[ C_d(\alpha_E)-C_d(\alpha_0)\right ] \Delta C_{d,f}'' \\
+       C_{l,\text{dyn}}(t) &= C_{l,\text{circ}} + \pi T_u \omega   \\
+   %   C_{d,\text{dyn}}(t) &=  C_d(\alpha_E) + (\alpha_{ac}-\alpha_E) C_{l,\text{dyn}} + \left[ C_d(\alpha_E)-C_d(\alpha_0)\right ] \Delta C_{d,f}'' \\
+      C_{d,\text{dyn}}(t)  &=  C_d(\alpha_E) + \left[(\alpha_{ac}-\alpha_E) +T_u \omega \right]C_{l,\text{circ}} + \left[ C_d(\alpha_E)-C_d(\alpha_0)\right ] \Delta C_{d,f}'' \\
    %     C_{m,\text{dyn}}(t) &=  C_m(\alpha_E) + C_{l,\text{dyn}} \Delta C_{m,f}'' - \frac{\pi}{2} T_u \omega\\
        C_{m,\text{dyn}}(t) &=  C_m(\alpha_E) - \frac{\pi}{2} T_u \omega\\
    \end{aligned}
@@ -245,7 +252,8 @@ with:
    \begin{aligned}
        \Delta C_{d,f}'' &= \frac{\sqrt{f_s^{st}(\alpha_E)}-\sqrt{x_4}}{2} - \frac{f_s^{st}(\alpha_E)-x_4}{4} 
    ,\qquad
-       x_4\ge 0
+       x_4\ge 0  \\
+    C_{l,\text{circ}}&= x_4 (\alpha_E-\alpha_0) C_{l,\alpha} +  (1-x_4) C_{l,{\text{fs}}}(\alpha_E) 
    \end{aligned}
 
 
@@ -258,7 +266,7 @@ Beddoes-Leishman 5-states model (UAMod=5)
 The 5-states (incompressible) dynamic stall model is similar to the Beddoes-Leishman 4-states model (UAMod=4), but 
 adds a 5th state to represent vortex generation. 
 It is enabled using ``UAMod=5``. The model uses :math:`C_n` and :math:`C_c` as main physical quantities.  
-Linearization of the model will be available in the future.
+Linearization of the model is available.
 
 
 
@@ -272,7 +280,7 @@ Oye model (UAMod=6)
 
 Oye's dynamic stall model is a one-state (continuous) model, formulated in :cite:`ad-Oye:1991` and described e.g. in :cite:`ad-Branlard:book`.
 The model attempts to capture trailing edge stall. 
-Linearization of the model will be available in the future.
+Linearization of the model is available.
 
 
 **State equation:**