diff --git a/docs/source/user/aerodyn/bibliography.bib b/docs/source/user/aerodyn/bibliography.bib index 248a77a523..cd9265b2dc 100644 --- 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 index f99fc1216c..49071ea999 100644 --- 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 index dfa632f1c9..16d095bc91 100644 --- 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 index 7c66917bb8..da05583613 100644 --- 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:** diff --git a/modules/aerodyn/src/AeroDyn.f90 b/modules/aerodyn/src/AeroDyn.f90 index 355595e88d..608d001603 100644 --- a/modules/aerodyn/src/AeroDyn.f90 +++ b/modules/aerodyn/src/AeroDyn.f90 @@ -1850,12 +1850,6 @@ subroutine SetInputsForBEMT(p, u, m, indx, errStat, errMsg) m%BEMT_u(indx)%Vx(j,k) = dot_product( tmp, x_hat ) ! normal component (normal to the plane, not chord) of the inflow velocity of the jth node in the kth blade m%BEMT_u(indx)%Vy(j,k) = dot_product( tmp, y_hat ) ! tangential component (tangential to the plane, not chord) of the inflow velocity of the jth node in the kth blade - - ! inputs for DBEMT (DBEMT_Mod == DBEMT_cont_tauConst) - if (allocated(m%BEMT_u(indx)%Vx_elast_dot)) then - m%BEMT_u(indx)%Vx_elast_dot(j,k) = dot_product( u%BladeMotion(k)%TranslationAcc(:,j), x_hat ) ! normal component (normal to the plane, not chord) of the inflow velocity of the jth node in the kth blade - m%BEMT_u(indx)%Vy_elast_dot(j,k) = dot_product( u%BladeMotion(k)%TranslationAcc(:,j), y_hat ) ! tangential component (tangential to the plane, not chord) of the inflow velocity of the jth node in the kth blade - end if ! inputs for CUA (and CDBEMT): m%BEMT_u(indx)%omega_z(j,k) = dot_product( u%BladeMotion(k)%RotationVel( :,j), m%WithoutSweepPitchTwist(3,:,j,k) ) ! rotation of no-sweep-pitch coordinate system around z of the jth node in the kth blade @@ -2635,21 +2629,17 @@ SUBROUTINE ValidateInputData( InitInp, InputFileData, NumBl, ErrStat, ErrMsg ) !.................. if (InitInp%Linearize) then if (InputFileData%AFAeroMod /= AFAeroMod_Steady) then -!bjj: REMOVE when linearization has been tested - call SetErrStat( ErrID_Fatal, 'Steady blade airfoil aerodynamics must be used for linearization. Set AFAeroMod=1.', ErrStat, ErrMsg, RoutineName ) - !if (InputFileData%UAMod /= UA_HGM) then - ! call SetErrStat( ErrID_Fatal, 'When AFAeroMod=2, UAMod must be 4 for linearization. Set AFAeroMod=1 or UAMod=4.', ErrStat, ErrMsg, RoutineName ) - !end if + if (InputFileData%UAMod /= UA_HGM .and. InputFileData%UAMod /= UA_HGMV .and. InputFileData%UAMod /= UA_OYE) then + call SetErrStat( ErrID_Fatal, 'When AFAeroMod=2, UAMod must be 4, 5, or 6 for linearization. Set AFAeroMod=1, or, set UAMod=4, 5, or 6.', ErrStat, ErrMsg, RoutineName ) + end if end if if (InputFileData%WakeMod == WakeMod_FVW) then !bjj: note: among other things, WriteOutput values will not be calculated properly in AD Jacobians if FVW this is allowed call SetErrStat( ErrID_Fatal, 'FVW cannot currently be used for linearization. Set WakeMod=0 or WakeMod=1.', ErrStat, ErrMsg, RoutineName ) else if (InputFileData%WakeMod == WakeMod_DBEMT) then -!bjj: when linearization has been tested - call SetErrStat( ErrID_Fatal, 'DBEMT cannot currently be used for linearization. Set WakeMod=0 or WakeMod=1.', ErrStat, ErrMsg, RoutineName ) - !if (InputFileData%DBEMT_Mod /= DBEMT_cont_tauConst) then - ! call SetErrStat( ErrID_Fatal, 'DBEMT requires the continuous formulation with constant tau1 for linearization. Set DBEMT_Mod=3 or set WakeMod to 0 or 1.', ErrStat, ErrMsg, RoutineName ) - !end if + if (InputFileData%DBEMT_Mod /= DBEMT_cont_tauConst) then + call SetErrStat( ErrID_Fatal, 'DBEMT requires the continuous formulation with constant tau1 for linearization. Set DBEMT_Mod=3 or set WakeMod to 0 or 1.', ErrStat, ErrMsg, RoutineName ) + end if end if end if @@ -4979,8 +4969,8 @@ SUBROUTINE RotGetOP( t, u, p, p_AD, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, do j=1,p%NumBlades ! size(x%BEMT%DBEMT%element,2) do i=1,p%NumBlNds ! size(x%BEMT%DBEMT%element,1) - do k=1,size(x%BEMT%DBEMT%element(i,j)%vind_dot) - x_op(index) = x%BEMT%DBEMT%element(i,j)%vind_dot(k) + do k=1,size(x%BEMT%DBEMT%element(i,j)%vind_1) + x_op(index) = x%BEMT%DBEMT%element(i,j)%vind_1(k) index = index + 1 end do end do @@ -4989,14 +4979,23 @@ SUBROUTINE RotGetOP( t, u, p, p_AD, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, end if if (p%BEMT%UA%lin_nx>0) then - do j=1,p%NumBlades ! size(x%BEMT%UA%element,2) - do i=1,p%NumBlNds ! size(x%BEMT%UA%element,1) - do k=1,4 !size(x%BEMT%UA%element(i,j)%x) !linearize only first 4 states (5th is vortex) - x_op(index) = x%BEMT%UA%element(i,j)%x(k) + if (p%BEMT%UA%UAMod==UA_OYE) then + do j=1,p%NumBlades ! size(x%BEMT%UA%element,2) + do i=1,p%NumBlNds ! size(x%BEMT%UA%element,1) + x_op(index) = x%BEMT%UA%element(i,j)%x(4) index = index + 1 end do end do - end do + else + do j=1,p%NumBlades ! size(x%BEMT%UA%element,2) + do i=1,p%NumBlNds ! size(x%BEMT%UA%element,1) + do k=1,4 !size(x%BEMT%UA%element(i,j)%x) !linearize only first 4 states (5th is vortex) + x_op(index) = x%BEMT%UA%element(i,j)%x(k) + index = index + 1 + end do + end do + end do + endif end if @@ -5032,8 +5031,8 @@ SUBROUTINE RotGetOP( t, u, p, p_AD, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, do j=1,p%NumBlades ! size(dxdt%BEMT%DBEMT%element,2) do i=1,p%NumBlNds ! size(dxdt%BEMT%DBEMT%element,1) - do k=1,size(dxdt%BEMT%DBEMT%element(i,j)%vind_dot) - dx_op(index) = dxdt%BEMT%DBEMT%element(i,j)%vind_dot(k) + do k=1,size(dxdt%BEMT%DBEMT%element(i,j)%vind_1) + dx_op(index) = dxdt%BEMT%DBEMT%element(i,j)%vind_1(k) index = index + 1 end do end do @@ -5042,14 +5041,23 @@ SUBROUTINE RotGetOP( t, u, p, p_AD, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, end if if (p%BEMT%UA%lin_nx>0) then - do j=1,p%NumBlades ! size(dxdt%BEMT%UA%element,2) - do i=1,p%NumBlNds ! size(dxdt%BEMT%UA%element,1) - do k=1,4 !size(dxdt%BEMT%UA%element(i,j)%x) don't linearize 5th state - dx_op(index) = dxdt%BEMT%UA%element(i,j)%x(k) + if (p%BEMT%UA%UAMod==UA_OYE) then + do j=1,p%NumBlades ! size(dxdt%BEMT%UA%element,2) + do i=1,p%NumBlNds ! size(dxdt%BEMT%UA%element,1) + dx_op(index) = dxdt%BEMT%UA%element(i,j)%x(4) index = index + 1 end do end do - end do + else + do j=1,p%NumBlades ! size(dxdt%BEMT%UA%element,2) + do i=1,p%NumBlNds ! size(dxdt%BEMT%UA%element,1) + do k=1,4 !size(dxdt%BEMT%UA%element(i,j)%x) don't linearize 5th state + dx_op(index) = dxdt%BEMT%UA%element(i,j)%x(k) + index = index + 1 + end do + end do + end do + endif end if call AD_DestroyRotContinuousStateType( dxdt, ErrStat2, ErrMsg2) @@ -5552,15 +5560,17 @@ SUBROUTINE Init_Jacobian_x( p, InitOut, ErrStat, ErrMsg) do j=1,p%NumBlades ! size(x%BEMT%DBEMT%element,2) do i=1,p%NumBlNds ! size(x%BEMT%DBEMT%element,1) NodeTxt = 'blade '//trim(num2lstr(j))//', node '//trim(num2lstr(i)) + if (p%BEMT%UA%UAMod/=UA_OYE) then - InitOut%LinNames_x(k) = 'x1 '//trim(NodeTxt)//', rad' - k = k + 1 + InitOut%LinNames_x(k) = 'x1 '//trim(NodeTxt)//', rad' + k = k + 1 - InitOut%LinNames_x(k) = 'x2 '//trim(NodeTxt)//', rad' - k = k + 1 - - InitOut%LinNames_x(k) = 'x3 '//trim(NodeTxt)//', -' - k = k + 1 + InitOut%LinNames_x(k) = 'x2 '//trim(NodeTxt)//', rad' + k = k + 1 + + InitOut%LinNames_x(k) = 'x3 '//trim(NodeTxt)//', -' + k = k + 1 + endif InitOut%LinNames_x(k) = 'x4 '//trim(NodeTxt)//', -' p%dx(k) = 0.001 ! x4 is a number between 0 and 1, so we need this to be small @@ -5729,17 +5739,23 @@ SUBROUTINE Perturb_x( p, n, perturb_sign, x, dx ) if (n <= p%BEMT%DBEMT%lin_nx) then - if (n <= p%BEMT%DBEMT%lin_nx/2) then ! x_p%BEMT%DBEMT%element(i,j)%vind, else x_p%BEMT%DBEMT%element(i,j)%vind_dot + if (n <= p%BEMT%DBEMT%lin_nx/2) then ! x_p%BEMT%DBEMT%element(i,j)%vind, else x_p%BEMT%DBEMT%element(i,j)%vind_1 call GetStateIndices( n, size(x%BEMT%DBEMT%element,2), size(x%BEMT%DBEMT%element,1), size(x%BEMT%DBEMT%element(1,1)%vind), Blade, BladeNode, StateIndex ) x%BEMT%DBEMT%element(BladeNode,Blade)%vind(StateIndex) = x%BEMT%DBEMT%element(BladeNode,Blade)%vind(StateIndex) + dx * perturb_sign else - call GetStateIndices( n - p%BEMT%DBEMT%lin_nx/2, size(x%BEMT%DBEMT%element,2), size(x%BEMT%DBEMT%element,1), size(x%BEMT%DBEMT%element(1,1)%vind_dot), Blade, BladeNode, StateIndex ) - x%BEMT%DBEMT%element(BladeNode,Blade)%vind_dot(StateIndex) = x%BEMT%DBEMT%element(BladeNode,Blade)%vind_dot(StateIndex) + dx * perturb_sign + call GetStateIndices( n - p%BEMT%DBEMT%lin_nx/2, size(x%BEMT%DBEMT%element,2), size(x%BEMT%DBEMT%element,1), size(x%BEMT%DBEMT%element(1,1)%vind_1), Blade, BladeNode, StateIndex ) + x%BEMT%DBEMT%element(BladeNode,Blade)%vind_1(StateIndex) = x%BEMT%DBEMT%element(BladeNode,Blade)%vind_1(StateIndex) + dx * perturb_sign endif else !call GetStateIndices( n - p%BEMT%DBEMT%lin_nx, size(x%BEMT%UA%element,2), size(x%BEMT%UA%element,1), size(x%BEMT%UA%element(1,1)%x), Blade, BladeNode, StateIndex ) - call GetStateIndices( n - p%BEMT%DBEMT%lin_nx, size(x%BEMT%UA%element,2), size(x%BEMT%UA%element,1), 4, Blade, BladeNode, StateIndex ) + + if (p%BEMT%UA%UAMod==UA_OYE) then + call GetStateIndices( n - p%BEMT%DBEMT%lin_nx, size(x%BEMT%UA%element,2), size(x%BEMT%UA%element,1), 1, Blade, BladeNode, StateIndex ) + StateIndex=4 ! Always the 4th one + else + call GetStateIndices( n - p%BEMT%DBEMT%lin_nx, size(x%BEMT%UA%element,2), size(x%BEMT%UA%element,1), 4, Blade, BladeNode, StateIndex ) + endif x%BEMT%UA%element(BladeNode,Blade)%x(StateIndex) = x%BEMT%UA%element(BladeNode,Blade)%x(StateIndex) + dx * perturb_sign end if @@ -5834,8 +5850,8 @@ SUBROUTINE Compute_dX(p, x_p, x_m, delta_p, delta_m, dX) do j=1,size(x_p%BEMT%DBEMT%element,2) ! number of blades do i=1,size(x_p%BEMT%DBEMT%element,1) ! number of nodes per blade - dX(indx_first:indx_first+1) = x_p%BEMT%DBEMT%element(i,j)%vind_dot - x_m%BEMT%DBEMT%element(i,j)%vind_dot - indx_first = indx_first + size(x_p%BEMT%DBEMT%element(i,j)%vind_dot) !+=2 + dX(indx_first:indx_first+1) = x_p%BEMT%DBEMT%element(i,j)%vind_1 - x_m%BEMT%DBEMT%element(i,j)%vind_1 + indx_first = indx_first + size(x_p%BEMT%DBEMT%element(i,j)%vind_1) !+=2 end do end do @@ -5843,12 +5859,21 @@ SUBROUTINE Compute_dX(p, x_p, x_m, delta_p, delta_m, dX) if (p%BEMT%UA%lin_nx>0) then - do j=1,size(x_p%BEMT%UA%element,2) ! number of blades - do i=1,size(x_p%BEMT%UA%element,1) ! number of nodes per blade - dX(indx_first:indx_first+3) = x_p%BEMT%UA%element(i,j)%x(1:4) - x_m%BEMT%UA%element(i,j)%x(1:4) - indx_first = indx_first + 4 ! = index_first += 4 + if (p%BEMT%UA%UAMod==UA_OYE) then + do j=1,size(x_p%BEMT%UA%element,2) ! number of blades + do i=1,size(x_p%BEMT%UA%element,1) ! number of nodes per blade + dX(indx_first) = x_p%BEMT%UA%element(i,j)%x(4) - x_m%BEMT%UA%element(i,j)%x(4) + indx_first = indx_first + 1 ! = index_first += 4 + end do end do - end do + else + do j=1,size(x_p%BEMT%UA%element,2) ! number of blades + do i=1,size(x_p%BEMT%UA%element,1) ! number of nodes per blade + dX(indx_first:indx_first+3) = x_p%BEMT%UA%element(i,j)%x(1:4) - x_m%BEMT%UA%element(i,j)%x(1:4) + indx_first = indx_first + 4 ! = index_first += 4 + end do + end do + endif end if diff --git a/modules/aerodyn/src/BEMT.f90 b/modules/aerodyn/src/BEMT.f90 index 4caee0b8f8..096e36ad0b 100644 --- a/modules/aerodyn/src/BEMT.f90 +++ b/modules/aerodyn/src/BEMT.f90 @@ -352,22 +352,6 @@ subroutine BEMT_AllocInput( u, p, errStat, errMsg ) end if u%Vy = 0.0_ReKi - if (p%DBEMT_Mod==DBEMT_cont_tauConst) then - allocate ( u%Vx_elast_dot( p%numBladeNodes, p%numBlades ), STAT = errStat2 ) - if ( errStat2 /= 0 ) then - call SetErrStat( ErrID_Fatal, 'Error allocating memory for u%Vx_dot.', errStat, errMsg, RoutineName ) - return - end if - u%Vx_elast_dot = 0.0_ReKi - - allocate ( u%Vy_elast_dot( p%numBladeNodes, p%numBlades ), STAT = errStat2 ) - if ( errStat2 /= 0 ) then - call SetErrStat( ErrID_Fatal, 'Error allocating memory for u%Vy_dot.', errStat, errMsg, RoutineName ) - return - end if - u%Vy_elast_dot = 0.0_ReKi - end if - allocate ( u%omega_z( p%numBladeNodes, p%numBlades ), STAT = errStat2 ) if ( errStat2 /= 0 ) then call SetErrStat( ErrID_Fatal, 'Error allocating memory for u%omega_z.', errStat, errMsg, RoutineName ) @@ -936,15 +920,6 @@ subroutine SetInputs_For_DBEMT(u_DBEMT, u, p, axInduction, tanInduction, Rtip) end do end do - if( allocated(u%Vx_elast_dot)) then ! only for DBEMT_Mod=DBEMT_cont_tauConst - do j = 1,p%numBlades - do i = 1,p%numBladeNodes - u_DBEMT%element(i,j)%vind_s_dot(1) = axInduction( i,j)*u%Vx_elast_dot(i,j) - u%omega_z(i,j)*tanInduction(i,j)*u%Vy(i,j) ! Eq. 41 - u_DBEMT%element(i,j)%vind_s_dot(2) = -tanInduction(i,j)*u%Vy_elast_dot(i,j) - u%omega_z(i,j)*axInduction( i,j)*u%Vx(i,j) ! Eq. 41 - end do - end do - end if - end subroutine SetInputs_For_DBEMT !.................................................................................................................................. @@ -2286,7 +2261,6 @@ subroutine WriteDEBUGValuesToFile(t, u, p, x, xd, z, OtherState, m, AFInfo) , "omega_z" & , "rLocal" , "UserProp" & , "AxInd", "TanInd" -! , "Vx_elast_dot" , "Vy_elast_dot" & end if @@ -2312,9 +2286,6 @@ subroutine WriteDEBUGValuesToFile(t, u, p, x, xd, z, OtherState, m, AFInfo) , u%UserProp( DEBUG_BLADENODE,DEBUG_BLADE) & , m%axInduction( DEBUG_BLADENODE,DEBUG_BLADE) & , m%tanInduction(DEBUG_BLADENODE,DEBUG_BLADE) -! these are not always allocated -! , u%Vx_elast_dot(DEBUG_BLADENODE,DEBUG_BLADE) & -! , u%Vy_elast_dot(DEBUG_BLADENODE,DEBUG_BLADE) & ! now write the residual function to a separate file: if ((DEBUG_nStep >= 0).AND.(DEBUG_nStep <= 450000).AND.(MOD(DEBUG_nStep,25) == 0)) then diff --git a/modules/aerodyn/src/BEMT_Registry.txt b/modules/aerodyn/src/BEMT_Registry.txt index e3b28be8f5..5c68c3bca9 100644 --- a/modules/aerodyn/src/BEMT_Registry.txt +++ b/modules/aerodyn/src/BEMT_Registry.txt @@ -155,8 +155,6 @@ typedef ^ ^ ReKi typedef ^ ^ ReKi TSR - - - "Tip-speed ratio (to check if BEM should be turned off)" - typedef ^ ^ ReKi Vx {:}{:} - - "Local axial velocity at node" m/s typedef ^ ^ ReKi Vy {:}{:} - - "Local tangential velocity at node" m/s -typedef ^ ^ ReKi Vx_elast_dot {:}{:} - - "Local relative axial acceleration at node (for CDBEMT)" "m/s^2" -typedef ^ ^ ReKi Vy_elast_dot {:}{:} - - "Local relative tangential acceleration at node (for CDBEMT)" "m/s^2" typedef ^ ^ ReKi omega_z {:}{:} - - "rotation of no-sweep-pitch-twist coordinate system around z (for CDBEMT and CUA)" "rad/s" typedef ^ ^ ReKi rLocal {:}{:} - - "Radial distance from center-of-rotation to node" m typedef ^ InputType ReKi Un_disk - - - "disk-averaged velocity normal to the rotor disk (for input to DBEMT)" m/s diff --git a/modules/aerodyn/src/BEMT_Types.f90 b/modules/aerodyn/src/BEMT_Types.f90 index bb5752eff9..0289ba8ffb 100644 --- a/modules/aerodyn/src/BEMT_Types.f90 +++ b/modules/aerodyn/src/BEMT_Types.f90 @@ -161,8 +161,6 @@ MODULE BEMT_Types REAL(ReKi) :: TSR !< Tip-speed ratio (to check if BEM should be turned off) [-] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: Vx !< Local axial velocity at node [m/s] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: Vy !< Local tangential velocity at node [m/s] - REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: Vx_elast_dot !< Local relative axial acceleration at node (for CDBEMT) [m/s^2] - REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: Vy_elast_dot !< Local relative tangential acceleration at node (for CDBEMT) [m/s^2] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: omega_z !< rotation of no-sweep-pitch-twist coordinate system around z (for CDBEMT and CUA) [rad/s] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: rLocal !< Radial distance from center-of-rotation to node [m] REAL(ReKi) :: Un_disk !< disk-averaged velocity normal to the rotor disk (for input to DBEMT) [m/s] @@ -4368,34 +4366,6 @@ SUBROUTINE BEMT_CopyInput( SrcInputData, DstInputData, CtrlCode, ErrStat, ErrMsg END IF DstInputData%Vy = SrcInputData%Vy ENDIF -IF (ALLOCATED(SrcInputData%Vx_elast_dot)) THEN - i1_l = LBOUND(SrcInputData%Vx_elast_dot,1) - i1_u = UBOUND(SrcInputData%Vx_elast_dot,1) - i2_l = LBOUND(SrcInputData%Vx_elast_dot,2) - i2_u = UBOUND(SrcInputData%Vx_elast_dot,2) - IF (.NOT. ALLOCATED(DstInputData%Vx_elast_dot)) THEN - ALLOCATE(DstInputData%Vx_elast_dot(i1_l:i1_u,i2_l:i2_u),STAT=ErrStat2) - IF (ErrStat2 /= 0) THEN - CALL SetErrStat(ErrID_Fatal, 'Error allocating DstInputData%Vx_elast_dot.', ErrStat, ErrMsg,RoutineName) - RETURN - END IF - END IF - DstInputData%Vx_elast_dot = SrcInputData%Vx_elast_dot -ENDIF -IF (ALLOCATED(SrcInputData%Vy_elast_dot)) THEN - i1_l = LBOUND(SrcInputData%Vy_elast_dot,1) - i1_u = UBOUND(SrcInputData%Vy_elast_dot,1) - i2_l = LBOUND(SrcInputData%Vy_elast_dot,2) - i2_u = UBOUND(SrcInputData%Vy_elast_dot,2) - IF (.NOT. ALLOCATED(DstInputData%Vy_elast_dot)) THEN - ALLOCATE(DstInputData%Vy_elast_dot(i1_l:i1_u,i2_l:i2_u),STAT=ErrStat2) - IF (ErrStat2 /= 0) THEN - CALL SetErrStat(ErrID_Fatal, 'Error allocating DstInputData%Vy_elast_dot.', ErrStat, ErrMsg,RoutineName) - RETURN - END IF - END IF - DstInputData%Vy_elast_dot = SrcInputData%Vy_elast_dot -ENDIF IF (ALLOCATED(SrcInputData%omega_z)) THEN i1_l = LBOUND(SrcInputData%omega_z,1) i1_u = UBOUND(SrcInputData%omega_z,1) @@ -4462,12 +4432,6 @@ SUBROUTINE BEMT_DestroyInput( InputData, ErrStat, ErrMsg ) IF (ALLOCATED(InputData%Vy)) THEN DEALLOCATE(InputData%Vy) ENDIF -IF (ALLOCATED(InputData%Vx_elast_dot)) THEN - DEALLOCATE(InputData%Vx_elast_dot) -ENDIF -IF (ALLOCATED(InputData%Vy_elast_dot)) THEN - DEALLOCATE(InputData%Vy_elast_dot) -ENDIF IF (ALLOCATED(InputData%omega_z)) THEN DEALLOCATE(InputData%omega_z) ENDIF @@ -4537,16 +4501,6 @@ SUBROUTINE BEMT_PackInput( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrMsg, Int_BufSz = Int_BufSz + 2*2 ! Vy upper/lower bounds for each dimension Re_BufSz = Re_BufSz + SIZE(InData%Vy) ! Vy END IF - Int_BufSz = Int_BufSz + 1 ! Vx_elast_dot allocated yes/no - IF ( ALLOCATED(InData%Vx_elast_dot) ) THEN - Int_BufSz = Int_BufSz + 2*2 ! Vx_elast_dot upper/lower bounds for each dimension - Re_BufSz = Re_BufSz + SIZE(InData%Vx_elast_dot) ! Vx_elast_dot - END IF - Int_BufSz = Int_BufSz + 1 ! Vy_elast_dot allocated yes/no - IF ( ALLOCATED(InData%Vy_elast_dot) ) THEN - Int_BufSz = Int_BufSz + 2*2 ! Vy_elast_dot upper/lower bounds for each dimension - Re_BufSz = Re_BufSz + SIZE(InData%Vy_elast_dot) ! Vy_elast_dot - END IF Int_BufSz = Int_BufSz + 1 ! omega_z allocated yes/no IF ( ALLOCATED(InData%omega_z) ) THEN Int_BufSz = Int_BufSz + 2*2 ! omega_z upper/lower bounds for each dimension @@ -4671,46 +4625,6 @@ SUBROUTINE BEMT_PackInput( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrMsg, END DO END DO END IF - IF ( .NOT. ALLOCATED(InData%Vx_elast_dot) ) THEN - IntKiBuf( Int_Xferred ) = 0 - Int_Xferred = Int_Xferred + 1 - ELSE - IntKiBuf( Int_Xferred ) = 1 - Int_Xferred = Int_Xferred + 1 - IntKiBuf( Int_Xferred ) = LBOUND(InData%Vx_elast_dot,1) - IntKiBuf( Int_Xferred + 1) = UBOUND(InData%Vx_elast_dot,1) - Int_Xferred = Int_Xferred + 2 - IntKiBuf( Int_Xferred ) = LBOUND(InData%Vx_elast_dot,2) - IntKiBuf( Int_Xferred + 1) = UBOUND(InData%Vx_elast_dot,2) - Int_Xferred = Int_Xferred + 2 - - DO i2 = LBOUND(InData%Vx_elast_dot,2), UBOUND(InData%Vx_elast_dot,2) - DO i1 = LBOUND(InData%Vx_elast_dot,1), UBOUND(InData%Vx_elast_dot,1) - ReKiBuf(Re_Xferred) = InData%Vx_elast_dot(i1,i2) - Re_Xferred = Re_Xferred + 1 - END DO - END DO - END IF - IF ( .NOT. ALLOCATED(InData%Vy_elast_dot) ) THEN - IntKiBuf( Int_Xferred ) = 0 - Int_Xferred = Int_Xferred + 1 - ELSE - IntKiBuf( Int_Xferred ) = 1 - Int_Xferred = Int_Xferred + 1 - IntKiBuf( Int_Xferred ) = LBOUND(InData%Vy_elast_dot,1) - IntKiBuf( Int_Xferred + 1) = UBOUND(InData%Vy_elast_dot,1) - Int_Xferred = Int_Xferred + 2 - IntKiBuf( Int_Xferred ) = LBOUND(InData%Vy_elast_dot,2) - IntKiBuf( Int_Xferred + 1) = UBOUND(InData%Vy_elast_dot,2) - Int_Xferred = Int_Xferred + 2 - - DO i2 = LBOUND(InData%Vy_elast_dot,2), UBOUND(InData%Vy_elast_dot,2) - DO i1 = LBOUND(InData%Vy_elast_dot,1), UBOUND(InData%Vy_elast_dot,1) - ReKiBuf(Re_Xferred) = InData%Vy_elast_dot(i1,i2) - Re_Xferred = Re_Xferred + 1 - END DO - END DO - END IF IF ( .NOT. ALLOCATED(InData%omega_z) ) THEN IntKiBuf( Int_Xferred ) = 0 Int_Xferred = Int_Xferred + 1 @@ -4896,52 +4810,6 @@ SUBROUTINE BEMT_UnPackInput( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, ErrStat, ErrMs END DO END DO END IF - IF ( IntKiBuf( Int_Xferred ) == 0 ) THEN ! Vx_elast_dot not allocated - Int_Xferred = Int_Xferred + 1 - ELSE - Int_Xferred = Int_Xferred + 1 - i1_l = IntKiBuf( Int_Xferred ) - i1_u = IntKiBuf( Int_Xferred + 1) - Int_Xferred = Int_Xferred + 2 - i2_l = IntKiBuf( Int_Xferred ) - i2_u = IntKiBuf( Int_Xferred + 1) - Int_Xferred = Int_Xferred + 2 - IF (ALLOCATED(OutData%Vx_elast_dot)) DEALLOCATE(OutData%Vx_elast_dot) - ALLOCATE(OutData%Vx_elast_dot(i1_l:i1_u,i2_l:i2_u),STAT=ErrStat2) - IF (ErrStat2 /= 0) THEN - CALL SetErrStat(ErrID_Fatal, 'Error allocating OutData%Vx_elast_dot.', ErrStat, ErrMsg,RoutineName) - RETURN - END IF - DO i2 = LBOUND(OutData%Vx_elast_dot,2), UBOUND(OutData%Vx_elast_dot,2) - DO i1 = LBOUND(OutData%Vx_elast_dot,1), UBOUND(OutData%Vx_elast_dot,1) - OutData%Vx_elast_dot(i1,i2) = ReKiBuf(Re_Xferred) - Re_Xferred = Re_Xferred + 1 - END DO - END DO - END IF - IF ( IntKiBuf( Int_Xferred ) == 0 ) THEN ! Vy_elast_dot not allocated - Int_Xferred = Int_Xferred + 1 - ELSE - Int_Xferred = Int_Xferred + 1 - i1_l = IntKiBuf( Int_Xferred ) - i1_u = IntKiBuf( Int_Xferred + 1) - Int_Xferred = Int_Xferred + 2 - i2_l = IntKiBuf( Int_Xferred ) - i2_u = IntKiBuf( Int_Xferred + 1) - Int_Xferred = Int_Xferred + 2 - IF (ALLOCATED(OutData%Vy_elast_dot)) DEALLOCATE(OutData%Vy_elast_dot) - ALLOCATE(OutData%Vy_elast_dot(i1_l:i1_u,i2_l:i2_u),STAT=ErrStat2) - IF (ErrStat2 /= 0) THEN - CALL SetErrStat(ErrID_Fatal, 'Error allocating OutData%Vy_elast_dot.', ErrStat, ErrMsg,RoutineName) - RETURN - END IF - DO i2 = LBOUND(OutData%Vy_elast_dot,2), UBOUND(OutData%Vy_elast_dot,2) - DO i1 = LBOUND(OutData%Vy_elast_dot,1), UBOUND(OutData%Vy_elast_dot,1) - OutData%Vy_elast_dot(i1,i2) = ReKiBuf(Re_Xferred) - Re_Xferred = Re_Xferred + 1 - END DO - END DO - END IF IF ( IntKiBuf( Int_Xferred ) == 0 ) THEN ! omega_z not allocated Int_Xferred = Int_Xferred + 1 ELSE @@ -6116,22 +5984,6 @@ SUBROUTINE BEMT_Input_ExtrapInterp1(u1, u2, tin, u_out, tin_out, ErrStat, ErrMsg END DO END DO END IF ! check if allocated -IF (ALLOCATED(u_out%Vx_elast_dot) .AND. ALLOCATED(u1%Vx_elast_dot)) THEN - DO i2 = LBOUND(u_out%Vx_elast_dot,2),UBOUND(u_out%Vx_elast_dot,2) - DO i1 = LBOUND(u_out%Vx_elast_dot,1),UBOUND(u_out%Vx_elast_dot,1) - b = -(u1%Vx_elast_dot(i1,i2) - u2%Vx_elast_dot(i1,i2)) - u_out%Vx_elast_dot(i1,i2) = u1%Vx_elast_dot(i1,i2) + b * ScaleFactor - END DO - END DO -END IF ! check if allocated -IF (ALLOCATED(u_out%Vy_elast_dot) .AND. ALLOCATED(u1%Vy_elast_dot)) THEN - DO i2 = LBOUND(u_out%Vy_elast_dot,2),UBOUND(u_out%Vy_elast_dot,2) - DO i1 = LBOUND(u_out%Vy_elast_dot,1),UBOUND(u_out%Vy_elast_dot,1) - b = -(u1%Vy_elast_dot(i1,i2) - u2%Vy_elast_dot(i1,i2)) - u_out%Vy_elast_dot(i1,i2) = u1%Vy_elast_dot(i1,i2) + b * ScaleFactor - END DO - END DO -END IF ! check if allocated IF (ALLOCATED(u_out%omega_z) .AND. ALLOCATED(u1%omega_z)) THEN DO i2 = LBOUND(u_out%omega_z,2),UBOUND(u_out%omega_z,2) DO i1 = LBOUND(u_out%omega_z,1),UBOUND(u_out%omega_z,1) @@ -6260,24 +6112,6 @@ SUBROUTINE BEMT_Input_ExtrapInterp2(u1, u2, u3, tin, u_out, tin_out, ErrStat, Er END DO END DO END IF ! check if allocated -IF (ALLOCATED(u_out%Vx_elast_dot) .AND. ALLOCATED(u1%Vx_elast_dot)) THEN - DO i2 = LBOUND(u_out%Vx_elast_dot,2),UBOUND(u_out%Vx_elast_dot,2) - DO i1 = LBOUND(u_out%Vx_elast_dot,1),UBOUND(u_out%Vx_elast_dot,1) - b = (t(3)**2*(u1%Vx_elast_dot(i1,i2) - u2%Vx_elast_dot(i1,i2)) + t(2)**2*(-u1%Vx_elast_dot(i1,i2) + u3%Vx_elast_dot(i1,i2)))* scaleFactor - c = ( (t(2)-t(3))*u1%Vx_elast_dot(i1,i2) + t(3)*u2%Vx_elast_dot(i1,i2) - t(2)*u3%Vx_elast_dot(i1,i2) ) * scaleFactor - u_out%Vx_elast_dot(i1,i2) = u1%Vx_elast_dot(i1,i2) + b + c * t_out - END DO - END DO -END IF ! check if allocated -IF (ALLOCATED(u_out%Vy_elast_dot) .AND. ALLOCATED(u1%Vy_elast_dot)) THEN - DO i2 = LBOUND(u_out%Vy_elast_dot,2),UBOUND(u_out%Vy_elast_dot,2) - DO i1 = LBOUND(u_out%Vy_elast_dot,1),UBOUND(u_out%Vy_elast_dot,1) - b = (t(3)**2*(u1%Vy_elast_dot(i1,i2) - u2%Vy_elast_dot(i1,i2)) + t(2)**2*(-u1%Vy_elast_dot(i1,i2) + u3%Vy_elast_dot(i1,i2)))* scaleFactor - c = ( (t(2)-t(3))*u1%Vy_elast_dot(i1,i2) + t(3)*u2%Vy_elast_dot(i1,i2) - t(2)*u3%Vy_elast_dot(i1,i2) ) * scaleFactor - u_out%Vy_elast_dot(i1,i2) = u1%Vy_elast_dot(i1,i2) + b + c * t_out - END DO - END DO -END IF ! check if allocated IF (ALLOCATED(u_out%omega_z) .AND. ALLOCATED(u1%omega_z)) THEN DO i2 = LBOUND(u_out%omega_z,2),UBOUND(u_out%omega_z,2) DO i1 = LBOUND(u_out%omega_z,1),UBOUND(u_out%omega_z,1) diff --git a/modules/aerodyn/src/DBEMT.f90 b/modules/aerodyn/src/DBEMT.f90 index 4d3eae6a17..3c12e5c3b4 100644 --- a/modules/aerodyn/src/DBEMT.f90 +++ b/modules/aerodyn/src/DBEMT.f90 @@ -17,6 +17,14 @@ ! See the License for the specific language governing permissions and ! limitations under the License. !********************************************************************************************************************************** +! +! References: +! [1] E. Branlard, B. Jonkman, G.R. Pirrung, K. Dixon, J. Jonkman (2022) +! Dynamic inflow and unsteady aerodynamics models for modal and stability analyses in OpenFAST, +! Journal of Physics: Conference Series, doi:10.1088/1742-6596/2265/3/032044 +! [2] R. Damiani, J.Jonkman +! DBEMT Theory Rev. 3 +! Unpublished module DBEMT use NWTC_Library @@ -208,7 +216,7 @@ subroutine DBEMT_Init( InitInp, u, p, x, OtherState, m, Interval, InitOut, ErrSt end if end do - p%lin_nx = p%numNodes*p%numBlades*4 ! vind and vind_dot + p%lin_nx = p%numNodes*p%numBlades*4 ! vind and vind_1 else p%lin_nx = 0 end if @@ -241,9 +249,8 @@ subroutine DBEMT_ReInit( p, x, OtherState, m ) do j=1,size(x%element,2) do i=1,size(x%element,1) - x%element(i,j)%vind = 0.0_ReKi - x%element(i,j)%vind_dot = 0.0_ReKi - x%element(i,j)%vind_1 = 0.0_ReKi + x%element(i,j)%vind = 0.0_ReKi ! Dynamic induced velocities + x%element(i,j)%vind_1 = 0.0_ReKi ! Reduced induced velocities end do end do @@ -306,10 +313,10 @@ subroutine DBEMT_InitStates( i, j, u, p, x, OtherState ) x%element(i,j)%vind(2) = u%element(i,j)%vind_s(2) if (p%DBEMT_Mod == DBEMT_cont_tauConst) then - x%element(i,j)%vind_dot(1) = u%element(i,j)%vind_s_dot(1) - x%element(i,j)%vind_dot(2) = u%element(i,j)%vind_s_dot(2) + x%element(i,j)%vind_1(1) = (1._ReKi - p%k_0ye)*u%element(i,j)%vind_s(1) ! Reduced velocity. Eq. (6) from [1] + x%element(i,j)%vind_1(2) = (1._ReKi - p%k_0ye)*u%element(i,j)%vind_s(2) else - x%element(i,j)%vind_1(1) = u%element(i,j)%vind_s(1) + x%element(i,j)%vind_1(1) = u%element(i,j)%vind_s(1) ! Intermediate velocity x%element(i,j)%vind_1(2) = u%element(i,j)%vind_s(2) end if @@ -455,7 +462,7 @@ subroutine ComputeTau1(u, p, m, tau1, errStat, errMsg) temp = (1.0-1.3*AxInd_disk)*Un_disk - tau1 = 1.1*u%R_disk/temp ! Eqn. 1.2 (note that we've eliminated possibility of temp being 0) + tau1 = 1.1*u%R_disk/temp ! Eq. (1) from [1] (note that we've eliminated possibility of temp being 0) tau1 = min(tau1, 100.0_ReKi) ! put a limit on this time constant so it isn't unrealistically long (particularly at initialization) end if @@ -484,8 +491,8 @@ subroutine ComputeTau2(i, j, u, p, tau1, tau2, k_tau_out) spanRatio = u%spanRatio end if - k_tau = 0.39 - 0.26*spanRatio**2 ! Eqn. 1.23b - tau2 = k_tau*tau1 ! Eqn. 1.7 or Eqn 1.23a + k_tau = 0.39 - 0.26*spanRatio**2 + tau2 = k_tau*tau1 ! Eq. (1) from [1] if (present(k_tau_out) ) k_tau_out = k_tau @@ -552,7 +559,8 @@ SUBROUTINE DBEMT_CalcContStateDeriv( i, j, t, u, p, x, OtherState, m, dxdt, ErrS ! LOCAL variables CHARACTER(*), PARAMETER :: RoutineName = 'DBEMT_CalcContStateDeriv' - REAL(ReKi) :: tauConst + REAL(ReKi) :: tau1inv + REAL(ReKi) :: tau2inv REAL(ReKi) :: tau1 REAL(ReKi) :: tau2 @@ -567,20 +575,15 @@ SUBROUTINE DBEMT_CalcContStateDeriv( i, j, t, u, p, x, OtherState, m, dxdt, ErrS call SetErrStat(ErrID_Fatal,"Continuous state derivatives cannot be calculated unless DBEMT_Mod is 3.",ErrStat,ErrMsg,RoutineName) return end if - tau1 = p%tau1_const - !call ComputeTau1( u, p, m, tau1, errStat, errMsg) call ComputeTau2(i, j, u, p, tau1, tau2) - - ! Implement Equation 37 from E.Branlard 16-Dec-2019 doc: + tau1inv = 1.0_ReKi/(tau1) + tau2inv = 1.0_ReKi/(tau2) - dxdt%vind = x%vind_dot - - tauConst = -1.0_ReKi/(tau1 * tau2) - - dxdt%vind_dot = tauConst * ( x%vind(:) + (tau1 + tau2)*x%vind_dot(:) & - - u%vind_s(:) - p%k_0ye*tau1*u%vind_s_dot(:) ) - + ! State derivatives, Eq. (7) from [1] + dxdt%vind_1 = -tau1inv * x%vind_1(:) + (1 - p%k_0ye) * tau1inv * u%vind_s(:) + dxdt%vind = tau2inv * x%vind_1(:) - tau2inv * x%vind(:) + p%k_0ye * tau2inv * u%vind_s(:) + END SUBROUTINE DBEMT_CalcContStateDeriv !---------------------------------------------------------------------------------------------------------------------------------- !> This subroutine implements the fourth-order Runge-Kutta Method (RK4) for numerically integrating ordinary differential equations: @@ -652,11 +655,11 @@ SUBROUTINE DBEMT_RK4( i, j, t, n, u, utimes, p, x, OtherState, m, ErrStat, ErrMs IF ( ErrStat >= AbortErrLev ) RETURN - k1%vind = p%dt * k1%vind - k1%vind_dot = p%dt * k1%vind_dot + k1%vind = p%dt * k1%vind + k1%vind_1 = p%dt * k1%vind_1 - x_tmp%vind = x%element(i,j)%vind + 0.5 * k1%vind - x_tmp%vind_dot = x%element(i,j)%vind_dot + 0.5 * k1%vind_dot + x_tmp%vind = x%element(i,j)%vind + 0.5 * k1%vind + x_tmp%vind_1 = x%element(i,j)%vind_1 + 0.5 * k1%vind_1 ! interpolate u to find u_interp = u(t + dt/2) TPlusHalfDt = t+0.5_DbKi*p%dt @@ -667,20 +670,20 @@ SUBROUTINE DBEMT_RK4( i, j, t, n, u, utimes, p, x, OtherState, m, ErrStat, ErrMs ! find xdot at t + dt/2 CALL DBEMT_CalcContStateDeriv( i, j, TPlusHalfDt, u_interp, p, x_tmp, OtherState, m, k2, ErrStat2, ErrMsg2 ) - k2%vind = p%dt * k2%vind - k2%vind_dot = p%dt * k2%vind_dot + k2%vind = p%dt * k2%vind + k2%vind_1 = p%dt * k2%vind_1 - x_tmp%vind = x%element(i,j)%vind + 0.5 * k2%vind - x_tmp%vind_dot = x%element(i,j)%vind_dot + 0.5 * k2%vind_dot + x_tmp%vind = x%element(i,j)%vind + 0.5 * k2%vind + x_tmp%vind_1 = x%element(i,j)%vind_1 + 0.5 * k2%vind_1 ! find xdot at t + dt/2 (note x_tmp has changed) CALL DBEMT_CalcContStateDeriv( i, j, TPlusHalfDt, u_interp, p, x_tmp, OtherState, m, k3, ErrStat2, ErrMsg2 ) - k3%vind = p%dt * k3%vind - k3%vind_dot = p%dt * k3%vind_dot + k3%vind = p%dt * k3%vind + k3%vind_1 = p%dt * k3%vind_1 - x_tmp%vind = x%element(i,j)%vind + k3%vind - x_tmp%vind_dot = x%element(i,j)%vind_dot + k3%vind_dot + x_tmp%vind = x%element(i,j)%vind + k3%vind + x_tmp%vind_1 = x%element(i,j)%vind_1 + k3%vind_1 ! interpolate u to find u_interp = u(t + dt) TPlusDt = t + p%dt @@ -691,11 +694,11 @@ SUBROUTINE DBEMT_RK4( i, j, t, n, u, utimes, p, x, OtherState, m, ErrStat, ErrMs ! find xdot at t + dt CALL DBEMT_CalcContStateDeriv( i, j, TPlusDt, u_interp, p, x_tmp, OtherState, m, k4, ErrStat2, ErrMsg2 ) - k4%vind = p%dt * k4%vind - k4%vind_dot = p%dt * k4%vind_dot + k4%vind = p%dt * k4%vind + k4%vind_1 = p%dt * k4%vind_1 - x%element(i,j)%vind = x%element(i,j)%vind + ( k1%vind + 2. * k2%vind + 2. * k3%vind + k4%vind ) / 6. - x%element(i,j)%vind_dot = x%element(i,j)%vind_dot + ( k1%vind_dot + 2. * k2%vind_dot + 2. * k3%vind_dot + k4%vind_dot ) / 6. + x%element(i,j)%vind = x%element(i,j)%vind + ( k1%vind + 2. * k2%vind + 2. * k3%vind + k4%vind ) / 6. + x%element(i,j)%vind_1 = x%element(i,j)%vind_1 + ( k1%vind_1 + 2. * k2%vind_1 + 2. * k3%vind_1 + k4%vind_1 ) / 6. END SUBROUTINE DBEMT_RK4 !---------------------------------------------------------------------------------------------------------------------------------- @@ -779,11 +782,11 @@ SUBROUTINE DBEMT_AB4( i, j, t, n, u, utimes, p, x, OtherState, m, ErrStat, ErrMs else - x%element(i,j)%vind = x%element(i,j)%vind + p%DT/24. * ( 55.*OtherState%xdot(1)%element(i,j)%vind - 59.*OtherState%xdot(2)%element(i,j)%vind & + x%element(i,j)%vind = x%element(i,j)%vind + p%DT/24. * ( 55.*OtherState%xdot(1)%element(i,j)%vind - 59.*OtherState%xdot(2)%element(i,j)%vind & + 37.*OtherState%xdot(3)%element(i,j)%vind - 9.*OtherState%xdot(4)%element(i,j)%vind ) - x%element(i,j)%vind_dot = x%element(i,j)%vind_dot + p%DT/24. * ( 55.*OtherState%xdot(1)%element(i,j)%vind_dot - 59.*OtherState%xdot(2)%element(i,j)%vind_dot & - + 37.*OtherState%xdot(3)%element(i,j)%vind_dot - 9.*OtherState%xdot(4)%element(i,j)%vind_dot ) + x%element(i,j)%vind_1 = x%element(i,j)%vind_1 + p%DT/24. * ( 55.*OtherState%xdot(1)%element(i,j)%vind_1 - 59.*OtherState%xdot(2)%element(i,j)%vind_1 & + + 37.*OtherState%xdot(3)%element(i,j)%vind_1 - 9.*OtherState%xdot(4)%element(i,j)%vind_1 ) endif @@ -861,13 +864,13 @@ SUBROUTINE DBEMT_ABM4( i, j, t, n, u, utimes, p, x, OtherState, m, ErrStat, ErrM IF ( ErrStat >= AbortErrLev ) RETURN - x%element(i,j)%vind = x_in%vind + p%DT/24. * ( 9. * xdot_pred%vind + 19. * OtherState%xdot(1)%element(i,j)%vind & + x%element(i,j)%vind = x_in%vind + p%DT/24. * ( 9. * xdot_pred%vind + 19. * OtherState%xdot(1)%element(i,j)%vind & - 5. * OtherState%xdot(2)%element(i,j)%vind & + 1. * OtherState%xdot(3)%element(i,j)%vind ) - x%element(i,j)%vind_dot = x_in%vind_dot + p%DT/24. * ( 9. * xdot_pred%vind_dot + 19. * OtherState%xdot(1)%element(i,j)%vind_dot & - - 5. * OtherState%xdot(2)%element(i,j)%vind_dot & - + 1. * OtherState%xdot(3)%element(i,j)%vind_dot ) + x%element(i,j)%vind_1 = x_in%vind_1 + p%DT/24. * ( 9. * xdot_pred%vind_1 + 19. * OtherState%xdot(1)%element(i,j)%vind_1 & + - 5. * OtherState%xdot(2)%element(i,j)%vind_1 & + + 1. * OtherState%xdot(3)%element(i,j)%vind_1 ) endif END SUBROUTINE DBEMT_ABM4 @@ -921,4 +924,4 @@ subroutine DBEMT_End( u, p, x, OtherState, m, ErrStat, ErrMsg ) END SUBROUTINE DBEMT_End -end module DBEMT \ No newline at end of file +end module DBEMT diff --git a/modules/aerodyn/src/DBEMT_Registry.txt b/modules/aerodyn/src/DBEMT_Registry.txt index c5fa20f529..21a726a068 100644 --- a/modules/aerodyn/src/DBEMT_Registry.txt +++ b/modules/aerodyn/src/DBEMT_Registry.txt @@ -29,8 +29,7 @@ typedef ^ ^ ReKi rLocal { typedef ^ InitOutputType ProgDesc Ver - - - "This module's name, version, and date" - typedef ^ DBEMT_ElementContinuousStateType R8Ki vind {2} - - "The filtered induced velocity, [1,i,j] is the axial induced velocity (-Vx*a) at node i on blade j and [2,i,j] is the tantential induced velocity (Vy*a')" m/s -typedef ^ DBEMT_ElementContinuousStateType R8Ki vind_dot {2} - - "Time derivative of the filtered induced velocity, x%vind in CCSD" "m/s^2" -typedef ^ DBEMT_ElementContinuousStateType R8Ki vind_1 {2} - - "The filtered intermediate induced velocity" "m/s" +typedef ^ DBEMT_ElementContinuousStateType R8Ki vind_1 {2} - - "The filtered reduced or intermediate induced velocity" "m/s" # ..... States .................................................................................................................... # Define continuous (differentiable) states here: @@ -73,7 +72,6 @@ typedef ^ ParameterType IntKi DBEMT_Mod # ..... Inputs .................................................................................................................... typedef ^ DBEMT_ElementInputType ReKi vind_s {2} - - "The unfiltered induced velocity, [1] is the axial induced velocity (-Vx*a) and [2] is the tangential induced velocity (Vy*a') at node i on blade j. Note that the inputs are used only operated on at a particular node and blade, so we don't store all elements" "m/s" -typedef ^ DBEMT_ElementInputType ReKi vind_s_dot {2} - - "The first time derivative of the unfiltered induced velocity, u%vind_s" "m/s^2" typedef ^ DBEMT_ElementInputType ReKi spanRatio - - - "Normalized span location of blade node" - # Define inputs that are contained on the mesh here: # diff --git a/modules/aerodyn/src/DBEMT_Types.f90 b/modules/aerodyn/src/DBEMT_Types.f90 index 04f0fb3af0..d87e22b1a4 100644 --- a/modules/aerodyn/src/DBEMT_Types.f90 +++ b/modules/aerodyn/src/DBEMT_Types.f90 @@ -54,8 +54,7 @@ MODULE DBEMT_Types ! ========= DBEMT_ElementContinuousStateType ======= TYPE, PUBLIC :: DBEMT_ElementContinuousStateType REAL(R8Ki) , DIMENSION(1:2) :: vind !< The filtered induced velocity, [1,i,j] is the axial induced velocity (-Vx*a) at node i on blade j and [2,i,j] is the tantential induced velocity (Vy*a') [m/s] - REAL(R8Ki) , DIMENSION(1:2) :: vind_dot !< Time derivative of the filtered induced velocity, x%vind in CCSD [m/s^2] - REAL(R8Ki) , DIMENSION(1:2) :: vind_1 !< The filtered intermediate induced velocity [m/s] + REAL(R8Ki) , DIMENSION(1:2) :: vind_1 !< The filtered reduced or intermediate induced velocity [m/s] END TYPE DBEMT_ElementContinuousStateType ! ======================= ! ========= DBEMT_ContinuousStateType ======= @@ -102,7 +101,6 @@ MODULE DBEMT_Types ! ========= DBEMT_ElementInputType ======= TYPE, PUBLIC :: DBEMT_ElementInputType REAL(ReKi) , DIMENSION(1:2) :: vind_s !< The unfiltered induced velocity, [1] is the axial induced velocity (-Vx*a) and [2] is the tangential induced velocity (Vy*a') at node i on blade j. Note that the inputs are used only operated on at a particular node and blade, so we don't store all elements [m/s] - REAL(ReKi) , DIMENSION(1:2) :: vind_s_dot !< The first time derivative of the unfiltered induced velocity, u%vind_s [m/s^2] REAL(ReKi) :: spanRatio !< Normalized span location of blade node [-] END TYPE DBEMT_ElementInputType ! ======================= @@ -559,7 +557,6 @@ SUBROUTINE DBEMT_CopyElementContinuousStateType( SrcElementContinuousStateTypeDa ErrStat = ErrID_None ErrMsg = "" DstElementContinuousStateTypeData%vind = SrcElementContinuousStateTypeData%vind - DstElementContinuousStateTypeData%vind_dot = SrcElementContinuousStateTypeData%vind_dot DstElementContinuousStateTypeData%vind_1 = SrcElementContinuousStateTypeData%vind_1 END SUBROUTINE DBEMT_CopyElementContinuousStateType @@ -610,7 +607,6 @@ SUBROUTINE DBEMT_PackElementContinuousStateType( ReKiBuf, DbKiBuf, IntKiBuf, Ind Db_BufSz = 0 Int_BufSz = 0 Db_BufSz = Db_BufSz + SIZE(InData%vind) ! vind - Db_BufSz = Db_BufSz + SIZE(InData%vind_dot) ! vind_dot Db_BufSz = Db_BufSz + SIZE(InData%vind_1) ! vind_1 IF ( Re_BufSz .GT. 0 ) THEN ALLOCATE( ReKiBuf( Re_BufSz ), STAT=ErrStat2 ) @@ -643,10 +639,6 @@ SUBROUTINE DBEMT_PackElementContinuousStateType( ReKiBuf, DbKiBuf, IntKiBuf, Ind DbKiBuf(Db_Xferred) = InData%vind(i1) Db_Xferred = Db_Xferred + 1 END DO - DO i1 = LBOUND(InData%vind_dot,1), UBOUND(InData%vind_dot,1) - DbKiBuf(Db_Xferred) = InData%vind_dot(i1) - Db_Xferred = Db_Xferred + 1 - END DO DO i1 = LBOUND(InData%vind_1,1), UBOUND(InData%vind_1,1) DbKiBuf(Db_Xferred) = InData%vind_1(i1) Db_Xferred = Db_Xferred + 1 @@ -686,12 +678,6 @@ SUBROUTINE DBEMT_UnPackElementContinuousStateType( ReKiBuf, DbKiBuf, IntKiBuf, O OutData%vind(i1) = REAL(DbKiBuf(Db_Xferred), R8Ki) Db_Xferred = Db_Xferred + 1 END DO - i1_l = LBOUND(OutData%vind_dot,1) - i1_u = UBOUND(OutData%vind_dot,1) - DO i1 = LBOUND(OutData%vind_dot,1), UBOUND(OutData%vind_dot,1) - OutData%vind_dot(i1) = REAL(DbKiBuf(Db_Xferred), R8Ki) - Db_Xferred = Db_Xferred + 1 - END DO i1_l = LBOUND(OutData%vind_1,1) i1_u = UBOUND(OutData%vind_1,1) DO i1 = LBOUND(OutData%vind_1,1), UBOUND(OutData%vind_1,1) @@ -1972,7 +1958,6 @@ SUBROUTINE DBEMT_CopyElementInputType( SrcElementInputTypeData, DstElementInputT ErrStat = ErrID_None ErrMsg = "" DstElementInputTypeData%vind_s = SrcElementInputTypeData%vind_s - DstElementInputTypeData%vind_s_dot = SrcElementInputTypeData%vind_s_dot DstElementInputTypeData%spanRatio = SrcElementInputTypeData%spanRatio END SUBROUTINE DBEMT_CopyElementInputType @@ -2023,7 +2008,6 @@ SUBROUTINE DBEMT_PackElementInputType( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrSt Db_BufSz = 0 Int_BufSz = 0 Re_BufSz = Re_BufSz + SIZE(InData%vind_s) ! vind_s - Re_BufSz = Re_BufSz + SIZE(InData%vind_s_dot) ! vind_s_dot Re_BufSz = Re_BufSz + 1 ! spanRatio IF ( Re_BufSz .GT. 0 ) THEN ALLOCATE( ReKiBuf( Re_BufSz ), STAT=ErrStat2 ) @@ -2056,10 +2040,6 @@ SUBROUTINE DBEMT_PackElementInputType( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrSt ReKiBuf(Re_Xferred) = InData%vind_s(i1) Re_Xferred = Re_Xferred + 1 END DO - DO i1 = LBOUND(InData%vind_s_dot,1), UBOUND(InData%vind_s_dot,1) - ReKiBuf(Re_Xferred) = InData%vind_s_dot(i1) - Re_Xferred = Re_Xferred + 1 - END DO ReKiBuf(Re_Xferred) = InData%spanRatio Re_Xferred = Re_Xferred + 1 END SUBROUTINE DBEMT_PackElementInputType @@ -2097,12 +2077,6 @@ SUBROUTINE DBEMT_UnPackElementInputType( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, Er OutData%vind_s(i1) = ReKiBuf(Re_Xferred) Re_Xferred = Re_Xferred + 1 END DO - i1_l = LBOUND(OutData%vind_s_dot,1) - i1_u = UBOUND(OutData%vind_s_dot,1) - DO i1 = LBOUND(OutData%vind_s_dot,1), UBOUND(OutData%vind_s_dot,1) - OutData%vind_s_dot(i1) = ReKiBuf(Re_Xferred) - Re_Xferred = Re_Xferred + 1 - END DO OutData%spanRatio = ReKiBuf(Re_Xferred) Re_Xferred = Re_Xferred + 1 END SUBROUTINE DBEMT_UnPackElementInputType @@ -2710,10 +2684,6 @@ SUBROUTINE DBEMT_ElementInputType_ExtrapInterp1(u1, u2, tin, u_out, tin_out, Err b = -(u1%vind_s(i1) - u2%vind_s(i1)) u_out%vind_s(i1) = u1%vind_s(i1) + b * ScaleFactor END DO - DO i1 = LBOUND(u_out%vind_s_dot,1),UBOUND(u_out%vind_s_dot,1) - b = -(u1%vind_s_dot(i1) - u2%vind_s_dot(i1)) - u_out%vind_s_dot(i1) = u1%vind_s_dot(i1) + b * ScaleFactor - END DO b = -(u1%spanRatio - u2%spanRatio) u_out%spanRatio = u1%spanRatio + b * ScaleFactor END SUBROUTINE DBEMT_ElementInputType_ExtrapInterp1 @@ -2778,11 +2748,6 @@ SUBROUTINE DBEMT_ElementInputType_ExtrapInterp2(u1, u2, u3, tin, u_out, tin_out, c = ( (t(2)-t(3))*u1%vind_s(i1) + t(3)*u2%vind_s(i1) - t(2)*u3%vind_s(i1) ) * scaleFactor u_out%vind_s(i1) = u1%vind_s(i1) + b + c * t_out END DO - DO i1 = LBOUND(u_out%vind_s_dot,1),UBOUND(u_out%vind_s_dot,1) - b = (t(3)**2*(u1%vind_s_dot(i1) - u2%vind_s_dot(i1)) + t(2)**2*(-u1%vind_s_dot(i1) + u3%vind_s_dot(i1)))* scaleFactor - c = ( (t(2)-t(3))*u1%vind_s_dot(i1) + t(3)*u2%vind_s_dot(i1) - t(2)*u3%vind_s_dot(i1) ) * scaleFactor - u_out%vind_s_dot(i1) = u1%vind_s_dot(i1) + b + c * t_out - END DO b = (t(3)**2*(u1%spanRatio - u2%spanRatio) + t(2)**2*(-u1%spanRatio + u3%spanRatio))* scaleFactor c = ( (t(2)-t(3))*u1%spanRatio + t(3)*u2%spanRatio - t(2)*u3%spanRatio ) * scaleFactor u_out%spanRatio = u1%spanRatio + b + c * t_out @@ -2902,14 +2867,6 @@ SUBROUTINE DBEMT_Input_ExtrapInterp1(u1, u2, tin, u_out, tin_out, ErrStat, ErrMs ENDDO DO i02 = LBOUND(u_out%element,2),UBOUND(u_out%element,2) DO i01 = LBOUND(u_out%element,1),UBOUND(u_out%element,1) - DO i1 = LBOUND(u_out%element(i01,i02)%vind_s_dot,1),UBOUND(u_out%element(i01,i02)%vind_s_dot,1) - b = -(u1%element(i01,i02)%vind_s_dot(i1) - u2%element(i01,i02)%vind_s_dot(i1)) - u_out%element(i01,i02)%vind_s_dot(i1) = u1%element(i01,i02)%vind_s_dot(i1) + b * ScaleFactor - END DO - ENDDO - ENDDO - DO i02 = LBOUND(u_out%element,2),UBOUND(u_out%element,2) - DO i01 = LBOUND(u_out%element,1),UBOUND(u_out%element,1) b = -(u1%element(i01,i02)%spanRatio - u2%element(i01,i02)%spanRatio) u_out%element(i01,i02)%spanRatio = u1%element(i01,i02)%spanRatio + b * ScaleFactor ENDDO @@ -2995,15 +2952,6 @@ SUBROUTINE DBEMT_Input_ExtrapInterp2(u1, u2, u3, tin, u_out, tin_out, ErrStat, E ENDDO DO i02 = LBOUND(u_out%element,2),UBOUND(u_out%element,2) DO i01 = LBOUND(u_out%element,1),UBOUND(u_out%element,1) - DO i1 = LBOUND(u_out%element(i01,i02)%vind_s_dot,1),UBOUND(u_out%element(i01,i02)%vind_s_dot,1) - b = (t(3)**2*(u1%element(i01,i02)%vind_s_dot(i1) - u2%element(i01,i02)%vind_s_dot(i1)) + t(2)**2*(-u1%element(i01,i02)%vind_s_dot(i1) + u3%element(i01,i02)%vind_s_dot(i1)))* scaleFactor - c = ( (t(2)-t(3))*u1%element(i01,i02)%vind_s_dot(i1) + t(3)*u2%element(i01,i02)%vind_s_dot(i1) - t(2)*u3%element(i01,i02)%vind_s_dot(i1) ) * scaleFactor - u_out%element(i01,i02)%vind_s_dot(i1) = u1%element(i01,i02)%vind_s_dot(i1) + b + c * t_out - END DO - ENDDO - ENDDO - DO i02 = LBOUND(u_out%element,2),UBOUND(u_out%element,2) - DO i01 = LBOUND(u_out%element,1),UBOUND(u_out%element,1) b = (t(3)**2*(u1%element(i01,i02)%spanRatio - u2%element(i01,i02)%spanRatio) + t(2)**2*(-u1%element(i01,i02)%spanRatio + u3%element(i01,i02)%spanRatio))* scaleFactor c = ( (t(2)-t(3))*u1%element(i01,i02)%spanRatio + t(3)*u2%element(i01,i02)%spanRatio - t(2)*u3%element(i01,i02)%spanRatio ) * scaleFactor u_out%element(i01,i02)%spanRatio = u1%element(i01,i02)%spanRatio + b + c * t_out diff --git a/modules/aerodyn/src/UnsteadyAero.f90 b/modules/aerodyn/src/UnsteadyAero.f90 index d6ea351670..abd633b8f8 100644 --- a/modules/aerodyn/src/UnsteadyAero.f90 +++ b/modules/aerodyn/src/UnsteadyAero.f90 @@ -18,6 +18,23 @@ ! limitations under the License. ! !********************************************************************************************************************************** +! References: +! +! [40] E. Branlard, B. Jonkman, G.R. Pirrung, K. Dixon, J. Jonkman (2022) +! Dynamic inflow and unsteady aerodynamics models for modal and stability analyses in OpenFAST, +! Journal of Physics: Conference Series, doi:10.1088/1742-6596/2265/3/032044 +! +! [41] E. Branlard, J. Jonkman, B.Jonkman (2020) +! Development plan for the aerodynamic linearization in OpenFAST +! Unpublished +! +! [70] User Documentation / AeroDyn / Unsteady Aerodynamics / Boing-Vertol model +! https://openfast.readthedocs.io/ +! +! [other] R. Damiani and G. Hayman (2017) +! The Unsteady Aerodynamics Module for FAST 8 +! NOTE: equations for this reference are labeled as x.y [n] where n is the number of the equation when several equations are given. + module UnsteadyAero ! This module uses equations defined in the document "The Unsteady Aerodynamics Module for FAST 8" by Rick Damiani and Greg Hayman, 28-Feb-2017 @@ -728,8 +745,10 @@ subroutine UA_SetParameters( dt, InitInp, p, AFInfo, AFIndx, ErrStat, ErrMsg ) p%Flookup = InitInp%Flookup p%ShedEffect = InitInp%ShedEffect - if (p%UAMod==UA_HGM) then + if (p%UAMod==UA_HGM .or. p%UAMod==UA_HGMV) then p%lin_nx = p%numBlades*p%nNodesPerBlade*4 ! 4 continuous states per node per blade (5th state isn't currently linearizable) + else if (p%UAMod==UA_OYE) then + p%lin_nx = p%numBlades*p%nNodesPerBlade*1 ! continuous state per node per blade, but stored at position 4 else p%lin_nx = 0 end if @@ -1648,7 +1667,7 @@ end function Failed end subroutine UA_UpdateDiscOtherState_BV !============================================================================== -!> Calculate angle of attacks using Boeing-Vertol model +!> Calculate angle of attacks using Boeing-Vertol model, see [70] !! Drag effective angle of attack needs extra computation subroutine BV_getAlphas(i, j, u, p, xd, BL_p, tc, alpha_34, alphaE_L, alphaLag_D, adotnorm) integer, intent(in ) :: i !< node index within a blade @@ -1707,7 +1726,7 @@ subroutine BV_getAlphas(i, j, u, p, xd, BL_p, tc, alpha_34, alphaE_L, alphaLag_D alphaLag_D = alpha_34 - dalphaD*isgn ! NOTE: not effective alpha yet for drag end subroutine BV_getAlphas !============================================================================== -!> Calculate gamma for lift and drag based rel thickness. See CACTUS BV_DynStall.f95 +!> Calculate gamma for lift and drag based rel thickness. See CACTUS BV_DynStall.f95 subroutine BV_getGammas(tc, umach, gammaL, gammaD) real(ReKi), intent(in) :: tc !< Relative thickness of airfoil real(ReKi), intent(in) :: umach !< Mach number of Urel, = Urel*MinfMinf (freestrem Mach), 0 for incompressible @@ -1750,7 +1769,7 @@ real(ReKi) function BV_TransA(BL_p) BV_TransA = .5_ReKi*dalphaMax ! transition region for fairing lagged AOA in pure lag model end function BV_TransA !============================================================================== -!> Calculate deltas to negative and postivive stall angle +!> Calculate deltas to negative and postivive stall angle, see [70] subroutine BV_delNP(adotnorm, alpha, alphaLag_D, BL_p, activeD, delN, delP) real(ReKi), intent(in) :: adotnorm !< alphadot * Tu real(ReKi), intent(in) :: alpha !< alpha (3/4) @@ -1782,7 +1801,7 @@ subroutine BV_delNP(adotnorm, alpha, alphaLag_D, BL_p, activeD, delN, delP) end if end subroutine BV_delNP !============================================================================== -!> Calculate effective angle of attack for drag coefficient, based on lagged angle of attack +!> Calculate effective angle of attack for drag coefficient, based on lagged angle of attack, see [70] real(ReKi) function BV_alphaE_D(adotnorm, alpha, alphaLag_D, BL_p, activeD) real(ReKi), intent(in) :: adotnorm !< alphadot * Tu real(ReKi), intent(in) :: alpha !< alpha (3/4) @@ -1811,7 +1830,7 @@ real(ReKi) function BV_alphaE_D(adotnorm, alpha, alphaLag_D, BL_p, activeD) end if end function BV_alphaE_D !============================================================================== -!> Activate dynamic stall for lift or drag +!> Activate dynamic stall for lift or drag, see [70] subroutine BV_UpdateActiveStates(adotnorm, alpha, alphaLag_D, alphaE_L, BL_p, activeL, activeD) real(ReKi), intent(in) :: adotnorm !< alphadot * Tu real(ReKi), intent(in) :: alpha !< alpha (3/4) @@ -2414,7 +2433,7 @@ SUBROUTINE HGM_Steady( i, j, u, p, x, AFInfo, ErrStat, ErrMsg ) x%x(3) = AFI_interp%Cl ! Not used elseif (p%UAMod==UA_HGM) then - x%x(3) = BL_p%c_lalpha * (alphaE-BL_p%alpha0) + x%x(3) = BL_p%c_lalpha * (alphaE-BL_p%alpha0) ! Clp ! calculate x%x(4) = fs_aF = f_st(alphaF): !alphaF = x%x(3)/BL_p%c_lalpha + BL_p%alpha0 ! p. 13 @@ -2513,7 +2532,7 @@ subroutine UA_CalcContStateDeriv( i, j, t, u_in, p, x, OtherState, AFInfo, m, dx ! find alphaF where FullyAttached(alphaF) = x(3) if (p%UAMod == UA_HGM) then !note: BL_p%c_lalpha cannot be zero. UA is turned off at initialization if this occurs. - alphaF = x%x(3)/BL_p%c_lalpha + BL_p%alpha0 ! p. 13 + alphaF = x%x(3)/BL_p%c_lalpha + BL_p%alpha0 ! Eq. 15 [40] else if (p%UAMod == UA_OYE) then alphaF = alpha_34 @@ -2550,8 +2569,8 @@ subroutine UA_CalcContStateDeriv( i, j, t, u_in, p, x, OtherState, AFInfo, m, dx call AddOrSub2Pi(real(x%x(1),ReKi), alpha_34) ! make sure we use the same alpha_34 for both x1 and x2 equations. if (p%ShedEffect) then - dxdt%x(1) = -1.0_R8Ki / Tu * (BL_p%b1 + p%c(i,j) * U_dot/(2*u%u**2)) * x%x(1) + BL_p%b1 * BL_p%A1 / Tu * alpha_34 - dxdt%x(2) = -1.0_R8Ki / Tu * (BL_p%b2 + p%c(i,j) * U_dot/(2*u%u**2)) * x%x(2) + BL_p%b2 * BL_p%A2 / Tu * alpha_34 + dxdt%x(1) = -1.0_R8Ki / Tu * (BL_p%b1 + p%c(i,j) * U_dot/(2*u%u**2)) * x%x(1) + BL_p%b1 * BL_p%A1 / Tu * alpha_34 ! Eq. 8 [40] + dxdt%x(2) = -1.0_R8Ki / Tu * (BL_p%b2 + p%c(i,j) * U_dot/(2*u%u**2)) * x%x(2) + BL_p%b2 * BL_p%A2 / Tu * alpha_34 ! Eq. 9 [40] else dxdt%x(1) = 0.0_ReKi dxdt%x(2) = 0.0_ReKi @@ -2560,8 +2579,8 @@ subroutine UA_CalcContStateDeriv( i, j, t, u_in, p, x, OtherState, AFInfo, m, dx if (p%UAMod == UA_HGM) then call AddOrSub2Pi(BL_p%alpha0, alphaE) Clp = BL_p%c_lalpha * (alphaE - BL_p%alpha0) + pi * Tu * u%omega ! Eq. 13 - dxdt%x(3) = -1.0_R8Ki / BL_p%T_p * x%x(3) + 1.0_ReKi / BL_p%T_p * Clp - dxdt%x(4) = -1.0_R8Ki / BL_p%T_f0 * x4 + 1.0_ReKi / BL_p%T_f0 * AFI_AlphaF%f_st + dxdt%x(3) = -1.0_R8Ki / BL_p%T_p * x%x(3) + 1.0_ReKi / BL_p%T_p * Clp ! Eq. 10 [40] + dxdt%x(4) = -1.0_R8Ki / BL_p%T_f0 * x4 + 1.0_ReKi / BL_p%T_f0 * AFI_AlphaF%f_st ! Eq. 11 [40] dxdt%x(5) = 0.0_R8Ki elseif (p%UAMod == UA_OYE) then @@ -3008,6 +3027,7 @@ subroutine UA_CalcOutput( i, j, t, u_in, p, x, xd, OtherState, AFInfo, y, misc, real(ReKi) :: cn_circ, tau_vl, tV_ratio real(ReKi) :: delta_c_df_primeprime real(ReKi), parameter :: delta_c_mf_primeprime = 0.0_ReKi + real(ReKi) :: cl_circ, cd_tors TYPE(UA_ElementContinuousStateType) :: x_in ! Continuous states at t ! for BV real(ReKi) :: alphaE_L, alphaE_D ! effective angle of attack for lift and drag @@ -3138,18 +3158,20 @@ subroutine UA_CalcOutput( i, j, t, u_in, p, x, xd, OtherState, AFInfo, y, misc, call AddOrSub2Pi(BL_p%alpha0, alphaE) cl_fa = (alphaE - BL_p%alpha0) * BL_p%c_lalpha - delta_c_df_primeprime = 0.5_ReKi * (sqrt(fs_aE) - sqrt(x4)) - 0.25_ReKi * (fs_aE - x4) ! Eq. 81 + delta_c_df_primeprime = 0.5_ReKi * (sqrt(fs_aE) - sqrt(x4)) - 0.25_ReKi * (fs_aE - x4) ! Eq. 20 [40] ! bjj: do we need to check that u%alpha is between -pi and + pi? - y%Cl = x4 * cl_fa + (1.0_ReKi - x4) * cl_fs + pi * Tu * u%omega ! Eq. 78 + cl_circ = x4 * cl_fa + (1.0_ReKi - x4) * cl_fs ! Eq. 19 [40] + y%Cl = cl_circ + pi * Tu * u%omega ! Eq. 16 [40] call AddOrSub2Pi(u%alpha, alphaE) - y%Cd = AFI_interp%Cd + (u%alpha - alphaE) * y%Cl + (AFI_interp%Cd - BL_p%Cd0) * delta_c_df_primeprime ! Eq. 79 + cd_tors = cl_circ * Tu * u%omega + y%Cd = AFI_interp%Cd + (alpha_34 - alphaE) * cl_circ + (AFI_interp%Cd - BL_p%Cd0) * delta_c_df_primeprime + cd_tors ! Eq. 17 [40] if (AFInfo%ColCm == 0) then ! we don't have a cm column, so make everything 0 y%Cm = 0.0_ReKi else - y%Cm = AFI_interp%Cm + y%Cl * delta_c_mf_primeprime - piBy2 * Tu * u%omega ! Eq. 80 + y%Cm = AFI_interp%Cm + y%Cl * delta_c_mf_primeprime - piBy2 * Tu * u%omega ! Eq. 18 [40] end if y%Cn = y%Cl*cos(u%alpha) + y%Cd*sin(u%alpha) @@ -3170,10 +3192,10 @@ subroutine UA_CalcOutput( i, j, t, u_in, p, x, xd, OtherState, AFInfo, y, misc, tau_vl = tau_vl / Tu ! make this non-dimensional (to compare with T_VL) tV_ratio = min(1.5_ReKi, tau_vl/BL_p%T_VL) - delta_c_df_primeprime = 0.5_ReKi * (sqrt(fs_aE) - sqrt(x4)) - 0.25_ReKi * (fs_aE - x4) ! Eq. 81 + delta_c_df_primeprime = 0.5_ReKi * (sqrt(fs_aE) - sqrt(x4)) - 0.25_ReKi * (fs_aE - x4) ! Eq. 20 [40] call AddOrSub2Pi(u%alpha, alphaE) - y%Cd = AFI_interp%Cd + (u%alpha - alphaE) * y%Cn + (AFI_interp%Cd - BL_p%Cd0) * delta_c_df_primeprime ! Eq. 79 + y%Cd = AFI_interp%Cd + (u%alpha - alphaE) * y%Cn + (AFI_interp%Cd - BL_p%Cd0) * delta_c_df_primeprime ! Eq. 79 [41] if (AFInfo%ColCm == 0) then ! we don't have a cm column, so make everything 0 @@ -3453,7 +3475,7 @@ subroutine UA_CalcOutput( i, j, t, u_in, p, x, xd, OtherState, AFInfo, y, misc, contains !> Calc Outputs for Boieng-Vertol dynamic stall - !! See BV_DynStall.f95 of CACTUS, notations kept more or less consistent + !! See BV_DynStall.f95 of CACTUS, and [70], notations kept more or less consistent subroutine BV_CalcOutput() real(ReKi) :: alpha_50 real(ReKi) :: Cm25_stat diff --git a/modules/aerodyn/src/UnsteadyAero_Driver.f90 b/modules/aerodyn/src/UnsteadyAero_Driver.f90 index cd1591a122..30dc3c0480 100644 --- a/modules/aerodyn/src/UnsteadyAero_Driver.f90 +++ b/modules/aerodyn/src/UnsteadyAero_Driver.f90 @@ -289,18 +289,47 @@ subroutine setUAinputs(n,u,t,dt,dvrInitInp,timeArr,AOAarr,Uarr,OmegaArr) real(ReKi), intent(in) :: OmegaArr(:) integer :: indx real(ReKi) :: phase + real(ReKi) :: d_ref2AC + real(ReKi) :: alpha_ref + real(ReKi) :: U_ref + real(ReKi) :: v_ref(2) + real(ReKi) :: v_34(2) + logical, parameter :: OscillationAtMidChord=.true. ! for legacy, use false + logical, parameter :: VelocityAt34 =.true. ! for legacy, use false u%UserProp = 0 u%Re = dvrInitInp%Re if ( dvrInitInp%SimMod == 1 ) then + if (OscillationAtMidChord) then + d_ref2AC =-0.25_ReKi ! -0.25: oscillations at mid_chord + else + d_ref2AC = 0.0_ReKi ! 0: oscillations at AC + endif + U_ref = dvrInitInp%InflowVel ! m/s + t = (n-1)*dt phase = (n+dvrInitInp%Phase-1)*2*pi/dvrInitInp%StepsPerCycle - u%alpha = (dvrInitInp%Amplitude * sin(phase) + dvrInitInp%Mean)*D2R ! This needs to be in radians - ! u%omega = dvrInitInp%Amplitude * cos(phase) * dvrInitInp%Frequency * pi**2 / 90.0 ! This needs to be in radians derivative: d_alpha /d_t + alpha_ref = (dvrInitInp%Amplitude * sin(phase) + dvrInitInp%Mean)*D2R ! This needs to be in radians + v_ref(1) = sin(alpha_ref)*U_ref + v_ref(2) = cos(alpha_ref)*U_ref u%omega = dvrInitInp%Amplitude * cos(phase) * 2*pi/dvrInitInp%StepsPerCycle / dt * D2R ! This needs to be in radians derivative: d_alpha /d_t - - u%U = dvrInitInp%InflowVel ! m/s + + u%v_ac(1) = v_ref(1) + u%omega * d_ref2AC* dvrInitInp%Chord + u%v_ac(2) = v_ref(2) + + v_34(1) = u%v_ac(1) + u%omega * 0.5* dvrInitInp%Chord + v_34(2) = u%v_ac(2) + + + u%alpha = atan2(u%v_ac(1), u%v_ac(2) ) ! + if (VelocityAt34) then + u%U = sqrt(v_34(1)**2 + v_34(2)**2) ! Using U at 3/4 + else + u%U = sqrt(u%v_ac(1)**2 + u%v_ac(2)**2) ! Using U at 1/4 + endif + + else indx = min(n,size(timeArr)) indx = max(1, indx) ! use constant data at initialization @@ -315,9 +344,9 @@ subroutine setUAinputs(n,u,t,dt,dvrInitInp,timeArr,AOAarr,Uarr,OmegaArr) elseif (n < 1) then t = (n-1)*dt end if + u%v_ac(1) = sin(u%alpha)*u%U + u%v_ac(2) = cos(u%alpha)*u%U end if - u%v_ac(1) = sin(u%alpha)*u%U - u%v_ac(2) = cos(u%alpha)*u%U end subroutine setUAinputs !---------------------------------------------------------------------------------------------------- diff --git a/modules/openfast-library/src/FAST_Lin.f90 b/modules/openfast-library/src/FAST_Lin.f90 index b0301a27f1..4645376c0f 100644 --- a/modules/openfast-library/src/FAST_Lin.f90 +++ b/modules/openfast-library/src/FAST_Lin.f90 @@ -5768,8 +5768,8 @@ SUBROUTINE PerturbOP(t, iLinTime, iMode, p_FAST, y_FAST, ED, BD, SrvD, AD, IfW, do j=1,size(AD%x(STATE_CURR)%rotors(1)%BEMT%DBEMT%element,2) do i=1,size(AD%x(STATE_CURR)%rotors(1)%BEMT%DBEMT%element,1) - indx_last = indx + size(AD%x(STATE_CURR)%rotors(1)%BEMT%DBEMT%element(i,j)%vind_dot) - 1 - call GetStateAry(p_FAST, iMode, t, AD%x(STATE_CURR)%rotors(1)%BEMT%DBEMT%element(i,j)%vind_dot, y_FAST%Lin%Modules(ThisModule)%Instance(1)%op_x_eig_mag( indx : indx_last), & + indx_last = indx + size(AD%x(STATE_CURR)%rotors(1)%BEMT%DBEMT%element(i,j)%vind_1) - 1 + call GetStateAry(p_FAST, iMode, t, AD%x(STATE_CURR)%rotors(1)%BEMT%DBEMT%element(i,j)%vind_1, y_FAST%Lin%Modules(ThisModule)%Instance(1)%op_x_eig_mag( indx : indx_last), & y_FAST%Lin%Modules(ThisModule)%Instance(1)%op_x_eig_phase(indx : indx_last) ) indx = indx_last + 1 end do diff --git a/reg_tests/CTestList.cmake b/reg_tests/CTestList.cmake index 34bd3772c7..f7051b3f8c 100644 --- a/reg_tests/CTestList.cmake +++ b/reg_tests/CTestList.cmake @@ -291,6 +291,7 @@ ad_regression("ad_VerticalAxis_OLAF" "aerodyn;bem") ad_regression("ad_BAR_CombinedCases" "aerodyn;bem") # NOTE: doing BAR at the end to avoid copy errors ad_regression("ad_BAR_OLAF" "aerodyn;bem") ad_regression("ad_BAR_SineMotion" "aerodyn;bem") +ad_regression("ad_BAR_SineMotion_UA4_DBEMT3" "aerodyn;bem") ad_regression("ad_BAR_RNAMotion" "aerodyn;bem") # BeamDyn regression tests diff --git a/reg_tests/r-test b/reg_tests/r-test index 1db296f230..84616202ad 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit 1db296f23054ed2580f0694395e5baaae5f66f4e +Subproject commit 84616202adc19515e08b5c78f92448aa848e5f12