Implementation of external forces in BeamDyn

[LaurenceWETI]

Hello everyone

@jjonkman @bjonkman

In a previous issue #548 we discussed how to implement the new loads resulting from shifting a mass along the blade in OpenFAST/ElastoDyn. The new loads are added in terms of external torque applied at the hub (HubPtLoad). This implementation allows to simulate the dynamic response of the wind turbine to changes in the direction of the shifted mass in- and outboard. However, the impact of the resulting Coriolis forces on the blade elements are not considered in that implementation. Therefore, the next step is to implement these forces in the source code of OpenFAST/BeamDyn. I started to get familiar with the theory manual and the source code of BeamDyn, and I would ask some questions relating to the procedure that I would go through to implement the resulting Coriolis forces.

The resulting Coriolis forces will be implanted as a distributed force vector in a separate Subroutine in BeamDyn.f90. The new distributed force will look like this u%DistrLoad%Force = u%DistrLoad%Force + u%DistrLoad%Coriolis The Coriolis force will be add only when the mass moved along the blade, else it equals zero u%DistrLoad%Coriolis = 0.0_ReKi The subroutine used to add the Coriolis force will be called in every subroutine uses the u%DistrLoad%Force i.e. in subroutines BD_GenerateDynamicElementAcc(); BD_GenerateDynamicElementForce();BD_Static();BD_GenerateStaticElement();BD_GenerateStaticElementForce();``BD_GenerateDynamicElementGA2() and BD_InputGlobalLocal() The parameters needed to compute the Coriolis force will be imported from a Simulink model through the S-Function as introduced in issue #548 and this work [1].

• Is it possible to consider the resulting Coriolis force as an external distributed force?

• Is the presented procedure feasible to be implemented in the source code of OpenFAST/BeamDyn?

After an internally proof of the modified code, I would pull a request in GitHub to get the changes in the source code reviewed from the OpenFAST developer.

I would appreciate any help.

Best regards

Laurence

[jjonkman]

Hi @LaurenceWETI,

Overall the approach sounds OK. Just a few comments:

• To have the external force passed from Simulink to BeamDyn, you'd have to add the external force as an input passed through Simulink, through the S-Function, through the OpenFAST library, and into BeamDyn.
• Presumably you are not proposing to share your implementation and include this feature into the official version of OpenFAST (or are you?); if so, there is no need to issue a PR.
• You mention the Coriolis force, but there are certainly other forces imposed by a the moving mass on a spinning blade. Do you intend to capture each of these forces?

Best regards,

[LaurenceWETI]

Hello @jjonkman,

Thank you for the quick reply.

• To pass signals from Simulink through the S-Function to BeamDyn I would follow the same procedures that I used by the addition of external torque applied at the hub (HubPtLoad). Part of the modified source code was discussed with @bjonkman and you in #548.

• The main point of PR is to make the modified codes publicly available on one hand, and on the hand, is to check the validity of these modifications by the OpenFAST developer. I think the functionalities of this feature must be first approved, before making a decision if this feature can be included into the official version of OpenFAST or not?

• Other forces imposed by moving the mass are intrinsic to the kinematics and kinetics of BeamDyn, i.e., as the element masses change the element forces (▁F ̂ ) will be dynamically updated to the new element masses. The only one force that can’t be computed is the Coriolis force resulting from change in the location of the mass along the blade. This is due to the facts that an element mass doesn’t change its location along the beam in the state-of-the-arte load simulation beam tools. Do you agree?

Best regards Laurence

[andrew-platt]

Hi @jjonkman and @LaurenceWETI,

I wonder if it might make sense to look at modifying the blade mounted structural controls to include Coriolis forces as needed? Those currently allow for full control by the DLL and could be extended to control through Simulink. It might be relatively straightforward to include the Coriolis force directly in those equations if it doesn't automatically fall out of the formulation and coupling.

Right now there is an assumption that blade mounted StCs are identical masses between blades, though they can be individually controlled. So this would take a little reworking of the StC interfacing if the masses need to be different. Some modifications would also need to be made if we prescribe displacement position rather than force (controller has access to the displacement along the blade, so it can control position in a closed loop manner).

One advantage of modifying the StC is that the coupling and most of the controller interfacing is already there. This would then work with both BeamDyn and ElastoDyn, and would require no modifications to BeamDyn. I think the only modifications needed would be within ServoDyn (for Simulink interface since that isn't complete) and the StC submodule (flag to indicate it is blade mounted, a distance vector from hub rotation axis, and rotor speed -- all of which ServoDyn already has).

Regards, Andy

[jjonkman]

Hi @andrew-platt and @LaurenceWETI,

Thanks for bringing up the structural control (StC) capability of OpenFAST. I should have brought this up yesterday as a different approach. The Coriolis force should already be included in this implementation, and as Andy mentioned, the motion of the mass can be controlled.

Best regards,

[LaurenceWETI]

Hello @andrew-platt and @jjonkman,

Thank you for pointing out the additional approach of the StCs. May I ask you some general questions regarding these controllers?

• Did I understand this correctly, The StCs are able to change the blade structural properties during a given simulation? Is there any theoretical description or a manual to use these controllers?

• The blade structural data stored in the FileInfo_In are those data that are derived from the BeamDyn or ElastoDyn blade input file? For example what is the InputFileData%StC_X_M ?

• If the blade element masses need to be changed, could I do this individually for each blade using the StCs?

• The controller used to change the blade element masses are developed in Simulink. Therefore, the StCs should be able to receive some control variables through the S-Function. Previously, I passed these variables through the S-Function and the OpenFAST glue code to ElastoDyn. Is that what do you mean with extending the StCs to be controlled through Simulink?

• What do you mean with this modification in StC submodule “flag to indicate it is blade mounted”?

Best regards, Laurence

[andrew-platt]

The StCs are essentially tuned mass dampers with the ability to control stiffness and damping through a controller. These are implemented through the ServoDyn module and may be mounted to blade, nacelle, substructure, or tower. The StC itself does not know where it is mounted as that is handled through the glue code. It is worth double checking the equations to verify that the Coriolis term is already included -- if it wasn't, some additional info may need to be passed along with a flag indicating it is blade mounted and needs the term.

The values for the blade StC are all stored within its input file. Documentation for it is here: https://openfast.readthedocs.io/en/dev/source/user/servodyn-stc/StC_index.html. There is also an update that includes updated controls interfacing for the StCs (this will be merged in shortly): https://github.com/OpenFAST/openfast/pull/664.

At present, all blade StCs are assumed to be the same constant mass, but may be motion controlled individually. Right now there is no provision for changing the mass during the simulation (this would take some reworking of the equations, but wouldn't be too difficult to do). Adding a mass change to the Simulink interface would be relatively straightforward -- there are two channels set aside in each StC channel group that could be used for mass changing (one for current mass value, and the other for the requested mass delta for the current timestep).

Regards, Andy

[jjonkman]

I'll just add that the Coriolis force is already included in the StC module, as mentioned in the theory documentation: https://openfast.readthedocs.io/en/dev/source/user/servodyn-stc/StC_Theory.html.

Best regards,

[LaurenceWETI]

Hello @andrew-platt and @jjonkman,

Thank you for the quick feedback and for your help. I read the Theory Manual for the Tuned Mass Damper Module in OpenFAST. and I tried to understand the source code of of the StCs. I see that there are some challenges to implement the principle of the hydro-pneumatic flywheel system (FW) as a TMD:

1. The TMD principle is based on a point mass that is attached to the structure at a specific location. Whereas, the FW system deals with distributed masses along the structure. This means that the attached mass StC_M needs to be changed from a point mass to a distributed masses at specific nodes along the structur.
2. The weight and the position of the distributed masses need to be dynamically and individually controlled for each blade within a given simulation.
3. The control parameters to determine the position and the weight of the distributed masses need to be passed through the S-Function, which is not available till now for the StCs, if I understand it correctly?

Can you give me please a feedback if these modifications are feasible in the source code of the StCs?

Best regards,

[andrew-platt]

A few thoughts on this:

1. For a distributed mass, you would need to define a line mesh option for the StC. This would take some additional code in both ServoDyn and in the glue-code to map a Line2 mesh. Ideally then the StC would define either a Line2 or point mesh, and the mapping would change depending on which one is committed. You would need to rework some of the equations in the StC code for the distributed mass -- this might be simplest as a new type of StC. Are you planning to do some kind of balloon/bladder for the liquid, and if so, how does that shape change depending on volume of mass added or inertial/gravitational loading on it?
2. It is possible to control the position of each blade mounted StC separately (in the control channel, enter 1,2,3 to assign blades 1-3 to different channels). The mass is not currently controlled, so this would require some modification -- there are two channels available per StC channel group for this (no changes to the DLL interface).
3. Correct. A limited portion of the StC channel groups are passed through the Simulink interface (not in dev yet, see my branch for this). You will need to expand this to include the additional channels.

Regards, Andy

[LaurenceWETI]

Hello @andrew-platt,

Thank you for your support.

I'm not sure that I could follow all the changes that you mentioned above to get the effects that I need in the StCs. I'm planning to control the dynamic behavior of the WT using the hydro-pneumatic flywheel system. the concept of such system is based on the change of the rotor inertia by moving a fluid mass between two accumulators in blade root and in blade tip, more details about this system can you find in our publications on ResearchGate.

I would do first the approach via BD, since I already did some modifications to pass the control variables through the Simulink interface. As soon as I finish the source code modifications, I would upload these modifications here to get the modified cods reviewed, if possible?

Best regards, Laurence

[LaurenceWETI]

Hello everyone,

@jjonkman @andrew-platt @bjonkman

I have two questions regarding adding forces to the applied distributed load %DistrLoad%Force.

The forces I want to add are the Coriolis forces resulting from shifting a fluid mass along the blade, see discussion above. The Coriolis force acting on the moving masses of the fluid, m_fluid, can be described by the equation below: (F_cor ) ⃗=2∙m_fluid∙(ω_rot ) ⃗∙(v_fluid ) ⃗ Where,(ω_rot ) ⃗is the angular velocity of the blade element and (v_fluid ) ⃗is the transitional fluid velocity.

This equation is implanted in BD source code to compute the x, y and z components of the distributed Coriolis force p%DistrCoriolisForce at node number jand beam member ias follow: X component: p%DistrCoriolisForce(1,j,i)=(2∙m_fluid∙(m%qp%vvv(5,j,i) * R2D)∙v_fluid)/((p%QPtWeight(j)*p%Jacobian(j,i))) Y component: p%DistrCoriolisForce(2,j,i)=(2∙m_fluid∙(m%qp%vvv(4,j,i) * R2D)∙v_fluid)/((p%QPtWeight(j)*p%Jacobian(j,i))) Z component: p%DistrCoriolisForce(3,j,i)=(2∙m_fluid∙(m%qp%vvv(6,j,i) * R2D)∙v_fluid)/((p%QPtWeight(j)*p%Jacobian(j,i)))

Where, (m%qp%vvv(4,j,i), (m%qp%vvv(5,j,i)and (m%qp%vvv(6,j,i) are the x, y and z component of the angular velocity expressed in l coordinates in rad/s, respectively. (p%QPtWeight(j)*p%Jacobian(j,i)) is the element length atj and i. The computed distributed Coriolis force components are added to the applied distributed forces u%DistrLoad%Force. u%DistrLoad%Force(:,j) = u%DistrLoad%Force(:,j) + p%DistrCoriolisForce(:,j,i)

The 5MW_Land_BD_DLL_WTurb.fst test is run without Aerodyn and ServoDy in order to exclude the effect of the additional forces. The p%DistrCoriolisForce are computed as expected, but they are not added to the applied distributed forces u%DistrLoad%Force, see figure below.

I see that in FATSolver .f90/SUBROUTINE BD_InputSolve, the applied distributed forces and moments are set to zero when AeroDyn is deactivated. could be that the reason why the Coriolis forces are not added to u%DistrLoad%Force? I tried to disable this condition, but the simulation crushed as it started to add the Coriolis forces to the applied distributed loads.

1. Did I made any mistake by implanting the Coriolis forces in the source code?

2. Is it possible to add forces to the applied distributed forces with deactivated aerodynamics?

I would appreciate any help.

Best regards

[jjonkman]

Dear @LaurenceWETI,

Where did you add the following equation: u%DistrLoad%Force(:,j) = u%DistrLoad%Force(:,j) + p%DistrCoriolisForce(:,j,i) in the source code? The inputs are sent to BeamDyn from the glue code in both the BD_UpdateStates() and BD_CalcOutput() subroutines; you may need to add the Coriolis force in both places.

FYI: I see that you declared the Coriolis force as a parameter (p%), but it is not really a parameter because it depends on more than just time (it appears to depend on states).

Best regards,

[LaurenceWETI]

Dear @jjonkman ,

Thank you very much for your quick reply.

That equation is implemented in SUBROUTINE BD_ChangeBladeMass(), which is called in BD_UpdateStates(). I’ll call BD_ChangeBladeMass() also in BD_CalcOutput() and see if that will make any change. Which data type would you suggest to declare the Coriolis force? As an input at time t, u, or as a misc/optimization variables,m?

Best regards

[jjonkman]

Dear @LaurenceWETI,

Yes, I agree that you'd have to add the Coriolis force to the input in both the BD_UpdateStates() and BD_CalcOutput() subroutines.

If the Coriolis force were an input to BeamDyn, the force should be set in the glue code rather than in BeamDyn. I'm guessing the Coriolis force should be a local variable. If the array size is large, it may make sense to make it a misc/optimization variable (m%) so that it is not allocated every subroutine call.

Best regards,

[LaurenceWETI]

Dear @jjonkman,

As you suggested I made the Coriolis force a misc/optimization variable(m%), and I called BD_ChangeBladeMass(), where the additional masses and the Coriolis forces are calculated, in both the BD_UpdateStates() and BD_CalcOutput() subroutines. Now the simulation crushed after circa 60 s as the rotor speed start to oscillate strongly, see figure below:

As before there are no changes in the distributed applied forces recognized. I think that this condition in FATSolver .f90/SUBROUTINE BD_InputSolve must be changed, to enable additional forces to be added to the distributed force, even when AerodDyn is deactivated:

I would set the calculation of the additional masses and the Coriolis forces in the glue code rather than BeamDyn, but the point is that this calculation is based on the mass matrix in BeamDyn, which is not accessible from the glue code.

Best regards,

[jjonkman]

Dear @LaurenceWETI,

I'm not sure why your model is going numerically unstable, but I would not think that moving the calculation of summing the Coriolis force with the aerodynamic force from BeamDyn to the glue code would eliminate the instability if the final resulting input used by BeamDyn is the same.

Regarding the instability, is the Coriolis force you are adding to the distributed load input smooth over time?

Best regards,

[LaurenceWETI]

Dear @jjonkman,

Thank you for your reply.

This is how the Coriolis force looks over the time.

It depends on the angular velocity and the varied fluid mass. The transitional velocity of the fluid is constant.

Best regards,

[jjonkman]

Dear @LaurenceWETI,

Is this just the Coriolis force at one node? Presumably you have additional forces at other nodes?

I'm a bit concerned that the non-smoothness of the force would cause numerical issues with the BeamDyn solver.

Best regards,

[LaurenceWETI]

Dear @jjonkman,

This is the distributed Coriolis force only at node 40. In this simulation the Coriolis force is applied from node 40 to node 48. See figure below.

FYI: As I called BD_ChangeBladeMass() in BD_CalcOutput(), I had to change the intent of the parameter data type from only (IN) to (INOUT) in BD_CalcOutput() subroutine and all subroutines that are called in BD_CalcOutput(), e.g. BD_JacobianPInput(). This is because of changing the intent of the mass matrix p%Mass0_QP() from (IN)to (INOUT), in order to add the fluid mass to the mass density in the mass matrix. Below is an example how I changed the mass matrix in a loop through the selected nodes (40 t0 48)

Best regards,

[jjonkman]

Dear @LaurenceWETI,

You should not change the INTENT of a parameter, which in the OpenFAST modularization framework are meant to be variables that are fully defined at initialization, with possible dependence only on time. Rather than changing the mass parameters directly, it would be better to calculate your additional mass and add it to the existing mass parameter(s) wherever they are needed in the solve.

Regardless, I'm guessing that is not the cause the instability. Again, I would be concerned that the instability is caused by the non-smoothness of the Coriolis force you've added. Can you implement a smoother form of your added force to verify that? Perhaps apply a simple function (like a sine wave) rather than the Coriolis force calculation you are implementing to verify that way you are applying additional forces works as expected.

Best regards,

[LaurenceWETI]

Dear @jinkman,

Thank you very much for support.

Your concerns were as you expected! After I implanted the Coriolis force in a smoother from (sine wave) my model goes numerically stable. However, I still can’t see the Coriolis forces added to the distributed applied load signal (NDFl).

Regarding the change of the mass parameters, I’ll check all subroutines where the mass parameter is needed and do the same loop to add the fluid mass to each beam/blade node. Which effect would you expect when this is done? And which effect has the change of the parameter INTENTon the model?

Best regards,

[jjonkman]

Dear @LaurenceWETI,

If you are not seeing your Coriolis force included in the BeamDyn output file, you are probably not adding it to the input early enough within BD_CalcOutput(). I would suggest adding your Coriolis force to the inputs at the very beginning of BD_UpdateStates() and BD_CalcOutput() (especially considering that you want the Coriolis force to be an input, but inputs should be calculated in the glue code, and you'd prefer to calculate it within BeamDyn).

Variables in the framework are given specific INTENTs to prevent errors in programming. Again, a parameter is meant to be fully definable at initialization and should not be changed at each time step. I'm not sure this applies to the specific parameter you want to change, but one problem I could see is that if other parameters depend on the parameter you change, you'll miss the propagation to those other parameters.

Best regards,

[LaurenceWETI]

Dear @jjonkman,

Thank you very much for your support.

I found a gap in my code that also contributed to the instability. The IF statement that I used to add the Coriolis force is only activated in case of changing the blade mass. Otherwise the solver take the control over the m%DistrCoriolisForce(:) in the rest of the simulation time, which made it much more unsmoothed. Therefore I set the Coriolis force to zero in the ELSE statement of that IF. That did help in the instability issue. In the example below I added only a distributed Coriolis force of 100 (N/m) in in-plan direction (only y component) on nodes from 40 to 48 without changing the blade element masses. The simulation completed successfully, see figure below.

• I think that the impact of additional force on the rotor speed is very small, don’t you? I tried to set the amplitude of the Coriolis force to greater value than 100 (N/m), but unfortunately the simulation crushed whenever the amplitude is greater than 100 (N/m)!

• As you suggested I add the Coriolis force in the very beginning of BD_UpdateStates() and BD_CalcOutput(). However, I still can’t see the Coriolis force in the output of BeamDyn. See the distributed applied load signal (NDFl) in the figure below.

Is this due to the disabled AeroDyn and ServoDyn that the solver prevent any changes in the applied distributed force? Is that related to this condition in FAT_Solver_ .f90/SUBROUTINE BD_InputSolve()?

Best regards,

[jjonkman]

Dear @LaurenceWETI,

This code snippet from SUBROUTINE BD_InputSolve() simply sets the applied distributed load input to BeamDyn to zero when AeroDyn is disabled. In your change to BeamDyn, you add your Coriolis force to this load input regardless of whether the original input is zero or nonzero, correct? Can you share how you add your Coriolis force to this input and which line(s) of the BeamDyn source code this addition takes place?

Best regards,

[LaurenceWETI]

Dear @jjonkman ,

Thank you for your reply.

It was not very clear to me whether this code snapped from SUBROUTINE BD_InputSolve() sets the applied distributed load to zero before or after the Coriolis force addition in BD source code. Now I understand that this code sets the input to BD to zero, which is changed after that in BD source due to the addition of Coriolis force. Thank you for the clarification.

Below I shared snap shots form BD of how I add the Coriolis force and the line(s) of the BD source coder where this addition takes place. The line numbers next to the subroutine name below refer to the place in the origin BD code. This is due to the additions that I did in SUBROUTINE BD_Init(), subroutine SetParameters() and subroutine Init_MiscVars(), which leads to shift the lines in my code for the origin one:

• SUBROUTINE BD_UpdateStates(), Line 1920:

• SUBROUTINE BD_CalcOutput(), Line 1964:

Where u%ExternalDeltaChargeState is the input form Simulink, which determines when the mass changes and thus Coriolis force must be added. This input is unequal zero when the mass start to change. BD_ChangeBladeMass()is the subroutine where the Coriolis force, m%DistrCoriolisForce, is added to the applied distributed load, u%DistrLoad%Force.

Best regards,

[jjonkman]

Dear @LaurenceWETI,

Can you share the full contents of SUBROUTINE BD_ChangeBladeMass()?

Best regards,

[LaurenceWETI]

Dear @jjonkman

Thank you for your quick reply.

Below is the full contenst of SUBROUTINE BD_ChangeBladeMass(). the first section of calculation and addition of fluid mass is deactivated. I will add the masses as you suggested wherever they are needed in the solve, rather than changing the mass density directly. I also deactivated the calculation of the Coriolis force components based on the resulted fluid mass. This is to be able to control the amplitude of the Coriolis force easily from Simulinkt with the input variable u%AmplitudeCoriolisForce

!-----------------------------------------------------------------------------------------------------------------------------------
!> LaurenceWETI: This subroutine calcuates and addes fluid mass and Coriolis force to the mass density per span unit 'm' and the applied distributed load u%DistrLoad%Force, respectively.
! The calcualtion is based on  u%ExternalChargeState and u%ExternalDeltaChargeState from Simulink.
SUBROUTINE BD_ChangeBladeMass(t,p,u,m,ErrStat,ErrMsg)

   REAL(DbKi),                      INTENT(IN   )  :: t                 !< Current simulation time in seconds
   TYPE(BD_InputType),              INTENT(INOUT)  :: u                 !< Inputs at t
   TYPE(BD_ParameterType),          INTENT(INOUT)  :: p                 !< Parameters
   TYPE(BD_MiscVarType),            INTENT(INOUT)  :: m                 !< misc/optimization variables

   INTEGER(IntKi),                  INTENT(  OUT)  :: ErrStat           !< Error status of the operation
   CHARACTER(*),                    INTENT(  OUT)  :: ErrMsg            !< Error message if ErrStat /= ErrID_None

   ! local variables
   INTEGER(IntKi)              :: nelem                                 !< index to the element counter
   REAL(BDKi)                  :: TipAcc_length(p%elem_total)           !< cumulative length of the tip accumulator - computed based on FE quadrature
   INTEGER(IntKi)              :: J_start                               !< First blade element of the tip accumulator close to blade root
   INTEGER(IntKi)              :: J_end                                 !< Last blade element of the tip accumulator close to blade tip
   REAL(BDKi)                  :: pi_roh_D                              !< pi/4 * fluid density * (diameter of tip accumulator)^2
   INTEGER(IntKi)              :: i
   INTEGER(IntKi)              :: j
   INTEGER(IntKi)              :: ErrStat2                              ! Temporary Error status
   CHARACTER(ErrMsgLen)        :: ErrMsg2                               ! Temporary Error message
   REAL(BDKi),dimension (6,6)  :: elem_ONE = (/1,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/)
   INTEGER(IntKi)              :: idx_qp
   CHARACTER(*), PARAMETER     :: RoutineName = 'BD_ChangeBladeMass'

    J_start                 = 40.0_BDKi     ! First blade element of the tip accumulator 
    J_end                   = 48.0_BDKi     ! Last blade element of the tip accumulator
    pi_roh_D                = 200           ! user defined parameter (example)

   ! Caculation of fluid mass  
!   TipAcc_length(:)        = 0.00_BDKi     ! Initialization
!   DO i=1,p%elem_total !  
!       DO j=J_start,J_end ! loop over tip accumulator nodes 
!           TipAcc_length(i)            = TipAcc_length(i) + p%QPtWeight(j)*p%Jacobian(j,i)
!          p%Mass0_QP(1:6,1:6,j)    = p%Mass0_QP(1:6,1:6,j) - elem_ONE * (p%BElmntFluidMass(j,i)/(p%QPtWeight(j)*p%Jacobian(j,i)))
!
!          IF       (TipAcc_length(i) <= u%ExternalChargeState * p%TipAcc_totoal_length) THEN
!                   p%BElmntFluidMass(j,i)  = p%QPtWeight(j)*p%Jacobian(j,i) * pi_roh_D
!                   p%Mass0_QP(1:6,1:6,j)   = p%Mass0_QP(1:6,1:6,j) + elem_ONE * (p%BElmntFluidMass(j,i)/(p%QPtWeight(j)*p%Jacobian(j,i)))
!                   p%BElmntMass(j,i)       = p%Mass0_QP(1,1,j) * p%QPtWeight(j)*p%Jacobian(j,i)
!          ELSEIF   ((TipAcc_length(i) > u%ExternalChargeState * p%TipAcc_totoal_length) .AND. ((TipAcc_length(i) - u%ExternalChargeState * p%TipAcc_totoal_length) <= p%QPtWeight(j)*p%Jacobian(j,i))) THEN
!                   p%BElmntFluidMass(j,i)  = (p%QPtWeight(j)*p%Jacobian(j,i) - (TipAcc_length(i)- u%ExternalChargeState * p%TipAcc_totoal_length)) * pi_roh_D
!                   p%Mass0_QP(1:6,1:6,j)   = p%Mass0_QP(1:6,1:6,j) + elem_ONE * (p%BElmntFluidMass(j,i)/(p%QPtWeight(j)*p%Jacobian(j,i)))
!                   p%BElmntMass(j,i)       = p%Mass0_QP(1,1,j) * p%QPtWeight(j)*p%Jacobian(j,i)            
!          ELSE
!                   p%BElmntFluidMass(j,i)  =  0.00
!                   p%Mass0_QP(1:6,1:6,j)   = p%Mass0_QP(1:6,1:6,j) + elem_ONE * (p%BElmntFluidMass(j,i)/(p%QPtWeight(j)*p%Jacobian(j,i)))
!                   p%BElmntMass(j,i)       = p%Mass0_QP(1,1,j) * p%QPtWeight(j)*p%Jacobian(j,i)            
!          ENDIF        
!       ENDDO
!   ENDDO

   ! Caculation of Coriolis force 
    TipAcc_length(:)        = 0.00_BDKi     ! Initialization
   DO i=1,p%elem_total !  
       DO j=J_start,J_end ! loop over tip accumulator nodes 
           TipAcc_length(i)         = TipAcc_length(i) + p%QPtWeight(j)*p%Jacobian(j,i)

           IF       (TipAcc_length(i) <= u%ExternalChargeState * p%TipAcc_totoal_length) THEN
!                   m%DistrCoriolisForce(1,j,i) = (2 * p%BElmntFluidMass(j,i) * (m%qp%vvv(5,j,i) * R2D) * (u%ExternalDeltaChargeState/p%dt * p%TipAcc_totoal_length )) / (p%QPtWeight(j)*p%Jacobian(j,i))
!                   m%DistrCoriolisForce(2,j,i) = (2 * p%BElmntFluidMass(j,i) * (m%qp%vvv(4,j,i) * R2D) * (u%ExternalDeltaChargeState/p%dt * p%TipAcc_totoal_length )) / (p%QPtWeight(j)*p%Jacobian(j,i))
!                   m%DistrCoriolisForce(3,j,i) = (2 * p%BElmntFluidMass(j,i) * (m%qp%vvv(6,j,i) * R2D) * (u%ExternalDeltaChargeState/p%dt * p%TipAcc_totoal_length )) / (p%QPtWeight(j)*p%Jacobian(j,i))
                    m%DistrCoriolisForce(1,j,i) = 0.0D0
                    m%DistrCoriolisForce(2,j,i) = u%AmplitudeCoriolisForce * u%ExternalDeltaChargeState 
                    m%DistrCoriolisForce(3,j,i) = 0.0D0
                    u%DistrLoad%Force(:,j)      = u%DistrLoad%Force(:,j) + m%DistrCoriolisForce(:,j,i) 

           ELSEIF   ((TipAcc_length(i) > u%ExternalChargeState * p%TipAcc_totoal_length) .AND. ((TipAcc_length(i) - u%ExternalChargeState * p%TipAcc_totoal_length) <= p%QPtWeight(j)*p%Jacobian(j,i))) THEN
                    m%DistrCoriolisForce(1,j,i) = 0.0D0
                    m%DistrCoriolisForce(2,j,i) = u%AmplitudeCoriolisForce * u%ExternalDeltaChargeState 
                    m%DistrCoriolisForce(3,j,i) = 0.0D0
                    u%DistrLoad%Force(:,j)      = u%DistrLoad%Force(:,j) + m%DistrCoriolisForce(:,j,i) 
           ELSE
                    m%DistrCoriolisForce(1,j,i) = 0.0D0
                    m%DistrCoriolisForce(2,j,i) = 0.0D0
                    m%DistrCoriolisForce(3,j,i) = 0.0D0
                    u%DistrLoad%Force(:,j)      = u%DistrLoad%Force(:,j) + m%DistrCoriolisForce(:,j,i) 

           ENDIF

       ENDDO
   ENDDO

   DO nelem=1,p%elem_total
      DO idx_qp=1,p%nqp
         p%qp%mmm(idx_qp,nelem)     =  p%Mass0_QP(1,1,(nelem-1)*p%nqp+idx_qp)  ! mass
      ENDDO
   ENDDO

      ! compute blade mass, CG, and IN for summary file:
   CALL BD_ComputeBladeMassNew( p, ErrStat2, ErrMsg2 )  !computes p%blade_mass,p%blade_CG,p%blade_IN
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

   RETURN

END SUBROUTINE BD_ChangeBladeMass

!-----------------------------------------------------------------------------------------------------------------------------------

Best regards,

[jjonkman]

Dear @LaurenceWETI,

It looks like you are changing u%DistrLoad%Force at the very beginning of SUBROUTINE BD_CalcOutput(), so, I'm not sure why you cannot see the Coriolis force in the BeamDyn write outputs. I would suggest stepping through the code with the debugger to figure out why the write outputs are not being updated properly.

This shouldn't effect the write outputs, but one issue I see in SUBROUTINE BD_UpdateStates() is that you only add the Coriolis force to the first element of the input array--that is, you only change u(1), whereas u in SUBROUTINE BD_UpdateStates() is an array covering multiple time steps at utimes.

Best regards,

[LaurenceWETI]

Dear @jjonkman,

Thank you for your help and your tips,

I just cloned the original OpenFAST code and wrote a very simple example in BD source code to add loads to the applied distributed force in the way that we discussed before. In the example below I added a 1000 (N/m) between 2 (s) and 3 (s) to u%DistrLoad%Force in both BD_UpdateStates()and BD_CalcOutput() subroutines:

• SUBROUTINE BD_UpdateStates(), Line 1921:

! LaurenceWETI: Add added a 1000 (N/m) between 2 (s) 3 (s) to the applied distributed force.
      IF( t >= 2 .AND. t <= 3  ) THEN
        u(1)%DistrLoad%Force  = 1000.0_ReKi
       ELSE
        u(1)%DistrLoad%Force  = 0.0_ReKi
      ENDIF

• SUBROUTINE BD_CalcOutput(), Line 1964:

! LaurenceWETI: Add added a 1000 (kN/m) between 2 (s) 3 (s) to the applied distributed force.
      IF( t >= 2 .AND. t <= 3  ) THEN
        u%DistrLoad%Force   = 1000.00_ReKi
       ELSE
        u%DistrLoad%Force   = 0.00_ReKi
      ENDIF

The figure below shows that the distributed applied load components at blade 1 station 1 change between the given time from 2 to 3 s. However, the changes are not in the same amplitude as it is given in the example of 1000 N/m! This simulation was run with disabled AeroDyn and ServoDyn.

I will add the modifications that I did in the source code before step by step and see which one causes the problem.

I have tow questions:

1. you mentioned in your previous comment that I add the Coriolis force only to the first element of the input array u(1). I tried to write a loop through to cover multiple time step at utimes, but it didn’t work. Can give an example how to fix this issue?
2. the amplitude of DNFl seems to be increased every time step! I would expect that the amplitude should be constant at 1000 N/m between 2 and 3 s. Do you have any explanation for this?

Best regards,

[jjonkman]

Dear @LaurenceWETI,

Regarding (1), instead of checking t in BD_UpdateStates(), I would check utimes and update the corresponding element of u.

Regarding (2), DistrLoad%Force is a 3-component vector in global coordinates. it looks like you are adding 1000 N/m to each coordinate. However, the output you are looking at is in a local coordinate system, which follows the motion (orientation) of the blade node. So, while the total magnitude should be the same (Magnitude = SQRT( 3*1000^2 )), the actual local components of force will change with the node's orientation.

Best regards,

[LaurenceWETI]

Dear @jjonkman,

Thank you very much for pointing that out.

Regarding (1): I’m still not sure if I understand how to update the corresponding element of u as an input of utimes. Could you please check if I update DistrLoad%Force to utimes correctly?

! Add 1000 (N/m) between 2 (s) 3 (s) to the applied distributed force.
        IF( utimes(t) >= 2 .AND. utimes(t) <= 3  ) THEN
                u(t)%DistrLoad%Force  = 1000.0_ReKi
        ELSE
                u(t)%DistrLoad%Force  = 0.0_ReKi
        ENDIF

Regarding (2), since, the Coriolis force depends on the 3-component vector of the angular velocity m%qp%vvv(4:6,j,i) which expressed in local coordinates, the resulting distributed Coriolis force m%DistrCoriolisForce cannot be added to u%DistrLoad%Force without transforming u%DistrLoad%Force from global to local or m%DistrCoriolisForce form local to global coordinates. As the Coriolis force must be added at the very beginning in the input of BD_UpdateStates() and BD_CalcOutput() subroutines, it makes more sense to transform m%DistrCoriolisForce rather than u%DistrLoad%Force. To transform the resulting Coriolis force from local to global I would multiply with the p%GlbRot matrix as follow:

   ! transform distributed Coriolis forces
   DO i=1,p%elem_total !  
       DO j=1,u%DistrLoad%Nnodes! 
          m%DistrCoriolisForce(1:3,j,i)  = MATMUL(m%DistrCoriolisForce(:,j,i),p%GlbRot)        
       ENDDO
   ENDDO

   ! add distributed Coriolis forces to applied distorted load
   DO i=1,p%elem_total !  
       DO j=1,u%DistrLoad%Nnodes! 
          u%DistrLoad%Force(:,j)        = u%DistrLoad%Force(:,j) + m%DistrCoriolisForce(:,j,i)         
       ENDDO
   ENDDO

Do you think this would work?

[jjonkman]

Dear @LaurenceWETI,

Regarding (1), in BD_UpdateStates(), u and utimes are arrays of size 2 or 3 (depending on InterpOrder in the OpenFAST primary input file). So, you must loop through the size of utimes, checking the value of each element along the way, and changing the corresponding element of u.

Regarding (2), I agree that you must convert m%DistrCoriolisForce to the global coordinate system before summing it with u%DistrLoad%Force, but you have not applied the correct transformation. The local coordinate system at each node will depend on that node. This orientation is an output of BeamDyn (y%BldMotion%Orientation). While this variable is not defined when you need to use it, you can see how it is calculated using calls to SUBROUTINEs BD_CrvCompose and BD_CrvMatrixR inBeamDyn_Subs.f90/Set_BldMotion.

Best regards,

[LaurenceWETI]

Dear @jjonkman,

Thank you for your quick reply.

I’m a bit concerned if I could apply the calculation of y%BldMotion%Orientation when I need. Therefore, I would try to calculate the Coriolis forces in Simulink based on the global components of the angular velocity NRVg from the output data. This has the advantage that the resulting Coriolis forces will be expressed in global coordinates and thus they can be added to the applied distributed load BD%Input(1,k)%DistrLoad%Force in FAST solver without the need to transform the coordinates. But I’m a bit concerned that the value of MAXInitINPUTSmust be extremely extended. This means for the 5MW_Land_BD_DLL_WTurb example, where the blades are divided in 49 stations, the MAXInitINPUTShas to be extended to 441 additional input channels to include the 3 components of Coriolis force for each blade.

1. Is there any upper limit for MAXInitINPUTS?
2. Do you agree with the procedure of implementation the Coriolis force in Simulink?

Best regards,

[jjonkman]

Dear @LaurenceWETI,

I don't really see a problem with increasing the number of inputs in Simulink to a large number.

That said, you had said earlier that you wanted to calculate the Coriolis force within BeamDyn so that you could use the mass matrix directly; this would be more complicated in Simulink. Also, if you add inputs to Simulink, you'll also have to modify the OpenFAST glue code so that the Coriolis forces from Simulink are added to the applied load inputs to BeamDyn.

While you can't use y%BldMotion%Orientation directly based on when you need to use it in the BeamDyn source code, as I mentioned before, you can calculate the orientation matrix using data that you do have available when you need it.

Best regards,

[LaurenceWETI]

Dear @jjonkman,

Thank you for your help, I am very grateful for your support.

I was not sure if I transformed the Coriolis force form local to global coordinates correctly. Therefore, I suggested to implant the Coriolis force in Simulink instead of BD. In the code attached below, I called Set_BldMotion_NoAcc() in BD_ChangeBladeMass() in order to calculate y%BldMotion%Orientation . Then I transformed the force components from local, m%DistrCoriolisForce, to global, m%DistrCoriolisForceG,.

!! calculate the orientation y%BldMotion%Orientation at each node 
CALL Set_BldMotion_NoAcc(p, x, m, y)

         do beta=1,p%nqp

               !-------------------------
               ! Applied distributed forces at Node 1 expressed in l, in N/m
            temp_vec = MATMUL(y%BldMotion%Orientation(:,:,beta), m%DistrCoriolisForce( :,beta))
            m%DistrCoriolisForceG(1,beta) = temp_vec(1)
            m%DistrCoriolisForceG(2,beta) = temp_vec(2)
            m%DistrCoriolisForceG(3,beta) = temp_vec(3)

         end do ! nodes

The resulting global components of the Coriolis force compared to the local components are presented in the figure below.

Do you have any idea how to proof if the transformation is correctly done? If I multiply the resulting global components again with the global matrix, MATMUL(m%DistrCoriolisForceG(:,i),p%GlbRot), should I become the origin local components of the Coriolis force, m%DistrCoriolisForce(:,i),?

FYI: In order to call BD_ChangeBladeMass(y) in BD_UpdateStatesI had to add the y data type to the attribute of BD_UpdateStatesin BD and in FAST_solve.

SUBROUTINE BD_UpdateStates( t, n, u, utimes, p, x, xd, z, OtherState, **y,** m, ErrStat, ErrMsg )

   REAL(DbKi),                      INTENT(IN   ) :: t          !< Current simulation time in seconds
   INTEGER(IntKi),                  INTENT(IN   ) :: n          !< Current simulation time step n = 0,1,...
   TYPE(BD_InputType),              INTENT(INOUT) :: u(:)       !< Inputs at utimes
   REAL(DbKi),                      INTENT(IN   ) :: utimes(:)  !< Times associated with u(:), in seconds
! @LaurenceWETI: Change p from only IN to INOUT
   TYPE(BD_ParameterType),          INTENT(IN   ) :: p          !< Parameters
   TYPE(BD_ContinuousStateType),    INTENT(INOUT) :: x          !< Input: Continuous states at t;
                                                                !!   Output: Continuous states at t + dt
   TYPE(BD_DiscreteStateType),      INTENT(INOUT) :: xd         !< Input: Discrete states at t;
                                                                !!   Output: Discrete states at t + dt
   TYPE(BD_ConstraintStateType),    INTENT(INOUT) :: z          !< Input: Initial guess of constraint states at t + dt;
                                                                !!   Output: Constraint states at t + dt
   TYPE(BD_OtherStateType),         INTENT(INOUT) :: OtherState !< Other states: Other states at t;
! @LaurenceWETI: Add y data type to call BD_ChangeBladeMass 
   **TYPE(BD_OutputType),           INTENT(INOUT) :: y          !< Outputs computed at t (Input only so that mesh con-**
                                                                !!   Output: Other states at t + dt
   TYPE(BD_MiscVarType),            INTENT(INOUT) :: m          !< misc/optimization variables
   INTEGER(IntKi),                  INTENT(  OUT) :: ErrStat    !< Error status of the operation
   CHARACTER(*),                    INTENT(  OUT) :: ErrMsg     !< Error message if ErrStat /= ErrID_None

! @LaurenceWETI: Call Subroutine BD_ChangeBladeMass when a change in the charge state takes place.
   INTEGER(IntKi)                               :: i           ! generic loop counter
   DO i=1,size(utimes)
      IF( u(i)%ExternalDeltaChargeState /= 0 ) THEN
       CALL BD_ChangeBladeMass(p,x,m,u(i),y,ErrStat,ErrMsg)
       ELSE
       m%DistrCoriolisForce = 0.0_ReKi
       m%DistrCoriolisForceG = 0.0_ReKi
      ENDIF
   ENDDO

Best regards,

[jjonkman]

Dear @LaurenceWETI,

BeamDyn output y%BldMotion%Orientation is a DCM, such that premultiplication of a vector in global coordinates by the DCM will transform the vector to a local coordinate system. So, to convert a vector in a local coordinate system to global, you should premultiply the vector by the transpose of the DCM (the transponse is the inverse for a DCM). To verify the implementation is correct, I would suggest generating a force that in local coordinates is constant and verify that the resulting output in local coordinates remains constant.

BeamDyn parameter p%GlbRot is not nodal dependent, and so, should not be used in the calculation.

In general I would recommend against changing arguments of the primary SUBROUTINEs in the OpenFAST modularization framework (_Init, _UpdateStates, _CalcOutput, _End). These arguments are dictated by the modularization framework and determine how we've structured the OpenFAST glue code. What I was suggesting above is not to call SUBROUTINE Set_BldMotion_NoAcc(), but rather re-implement the key equations needed to compute the DCM where you need them.

Best regards,

[LaurenceWETI]

Dear @jjonkman,

Thank you for pointing that out.

If I understand the interpolation approach of BD to deal with the finite rotation correclty, I’ve to call BD_CrvComposesubroutine twice: first to calculate the equivalent composition and then to interpolate the relative rotation. Subsequently, BD_CrvMatrixR() subroutine is called to compute the rotation tensor based on the resulting parameter rrfrom BD_CrvComposesubroutines.

But I’m not quite sure how to set the values of the two Winer-Milenkovic input parameters ppand qqto find the resulting composed parameter rrin SUBROUTINE BD_CrvCompose( rr, pp, qq, flag) to use them in my calculation in subroutine BD_ChangeBladeMass()?

Best regards,

[jjonkman]

Dear @LaurenceWETI,

I would call BD_CrvCompose() and BD_CrvMatrixR() just as they are called within Set_BldMotion_NoAcc(), e.g., with pp = x%q(4:6,temp_id).

Best regards,

[LaurenceWETI]

Dear @jjonkman,

Thank you very much for your reply.

The DCM y%BldMotion%Orientation in BeamDyn_Subs.f90/Set_BldMotion_NoAcc() subroutine is created in two cases: (case1) BD_MESH_FE and (case2) BE_MESH_QP. Since, the resulting Coriolis forces are applied at QPs, I generated and stored the DCM based on the second case as follow:

!! calculate the blade motion orientation tensor at each node 
   REAL(BDKi)                  :: BladeMotionOrientation(3,3,p%nqp)
   REAL(BDKi)                  :: BladeMotionOrientationTRANSPOSE(3,3,p%nqp)

   DO i=1,p%elem_total
      DO j=1,p%nqp
            temp_id2 = (i-1)*p%nqp + j + p%qp_indx_offset            ! Index to a node within element i           

            CALL BD_CrvCompose( cc, x%q(4:6,temp_id2), p%uuN0(4:6,j,i), FLAG_R1R2 )
            CALL BD_CrvCompose( cc0, p%Glb_crv, cc, FLAG_R1R2 )

               ! Create the DCM from the rotation parameters
            CALL BD_CrvMatrixR(cc0,temp_R)  ! returns temp_R (the transpose of the DCM orientation matrix)

               ! Store the DCM for the j'th node of the i'th element (in FAST coordinate system)            
                BladeMotionOrientation(1:3,1:3,temp_id2) = TRANSPOSE(temp_R)
                BladeMotionOrientationTRANSPOSE(1:3,1:3,temp_id2) = TRANSPOSE(BladeMotionOrientation(1:3,1:3,temp_id2))             

               ! Applied distributed Coriolis forces at Node j  expressed in g, in N/m
                m%DistrCoriolisForceG( :,j ) = MATMUL(BladeMotionOrientationTRANSPOSE(:,:,j), m%DistrCoriolisForce( :,j))

      END DO
   END DO

Where, m%DistrCoriolisForce is the calculated Coriolis force in local coordinates and m%DistrCoriolisForceG is the resulting output in global coordinates. The figure below shows these tow variables for a constant value of the y component of +10000 (N/m) between 1-2 s and -10000 (N/m) between 3-4 s at nodes from 40 to 48.

You mentioned in a previous #1003 that to verify the transformation the resulting force in global coordinates should be constant if the origin force in local coordinates was also constant. In this example the force in both coordinates are identical!

Do the transformed force components in the global coordinates look as you expected?

Best regards,

[jjonkman]

Dear @LaurenceWETI,

It looks like you mixed up the cases a bit. For case 2 (p%BldMotionNodeLoc = BD_MESH_QP), the first call to BD_CrvCompose() should use m%qp%uuu rather than x%q, and m%qp%uuu should be set by calling BD_RotationInterpQP().

For your test case, your setting what variable to plus or minus 10000?

Best regards,

[LaurenceWETI]

Dear @jjonkman,

Thank you for your reply.

Could you explain pleas why and where BD_RotationalInterpQP() subroutine need to called in my example? As you mentioned in a previous post BD_CrvCompose() and BD_CrvMatrixR() subroutines need to called as they are called in SUBROUTINE Set_BldMotion_NoAcc in BeamDyn_Subs.f90. Both subroutines are called in SUBROUTINE BD_ChangeBladeMassas follow:

!! calculate the blade motion orientation tensor at each node 
   DO i=1,p%elem_total
      DO j=1,p%nqp

         CALL BD_CrvCompose(cc,m%qp%uuu(4:6,j,i),p%uu0(4:6,j,i),FLAG_R1R2)   ! cc = m%qp%uuu(4:6,j,i) composed with p%uu0(4:6,j,i)
         CALL BD_CrvMatrixR(cc,m%qp%RR0(:,:,j,i))   ! returns RR0 (the transpose of the DCM orientation matrix)

               ! Store the DCM for the j'th node of the i'th element (in FAST coordinate system)            
                BladeMotionOrientation(1:3,1:3,j) = TRANSPOSE(m%qp%RR0(:,:,j,i))
                BladeMotionOrientationTRANSPOSE(1:3,1:3,j) = TRANSPOSE(BladeMotionOrientation(1:3,1:3,j))               

               ! Applied distributed Coriolis forces at Node j  expressed in g, in N/m
                m%DistrCoriolisForceG( :,j ) = MATMUL(BladeMotionOrientationTRANSPOSE(:,:,j), m%DistrCoriolisForce( :,j))
!               u%DistrLoad%Force(1:3,j)        = u%DistrLoad%Force(1:3,j) + m%DistrCoriolisForceG( 1:3,j ) 

      END DO
   END DO

The resulting transformed forces based on my transformation are not constant.

Regarding your question, the amplitude of Coriolis force is set to 10000 in Simulink as an external signal. This amplitude is multiplied by the derivative of the charge state, see subplot 1 in the figure above. The derivative of this signal is +1 between 2 and 3 s and -1 between 4 and 5.

Best regards,

[jjonkman]

Dear @LaurenceWETI,

It looks like you've now fixed the use of m%qp%uuu rather than x%q, but now you've dropped the second call to BD_CrvCompose involving p%Glb_crv and you've replaced the temp_R with m%qp%RR0.

The call to BD_RotationInterpQP() is needed because that is the SUBROUTINE that is setting m%qp%uuu.

To be honest, I'm not familiar with the BeamDyn source code (although I do understand the theory basis a bit). I'm just trying to get you to take what is already implement in the BeamDyn source code so that you can apply that code to your case.

Best regards,

[LaurenceWETI]

Dear @jjonkman,

Thank you for your reply.

I found the problem why the write outputs are not being updated properly. When additional outputs are added in OutListParameters.xlsx, e.g. Coriolis forces in my case, the generated subroutine SetOutParam() by Write_ChckOutLst.m/BeamDyn_SetOutParam.f90 doesn’t contain the channels validation check that is used in the original BeamDy_IO.f90. Below is the ifcondition at line 1454 in BeamDy_IO.f90 that must be added when Write_ChckOutLst.m is executed for additional output chanle/s:

      IF ( Indx > 0 ) THEN ! we found the channel name
         p%OutParam(I)%Indx     = ParamIndxAry(Indx)
         IF ( InvalidOutput( ParamIndxAry(Indx) ) ) THEN  ! but, it isn't valid for these settings
            p%OutParam(I)%Units = "INVALID"
            p%OutParam(I)%SignM = 0
         ELSE
            p%OutParam(I)%Units = ParamUnitsAry(Indx) ! it's a valid output

            if ( p%OutParam(I)%Indx >= N1DFxl .and. p%OutParam(I)%Indx <= N9DMzl ) p%OutInputs = .true.
         END IF
      ELSE ! this channel isn't valid
         p%OutParam(I)%Indx  = Time                 ! pick any valid channel (I just picked "Time" here because it's universal)
         p%OutParam(I)%Units = "INVALID"
         p%OutParam(I)%SignM = 0                    ! multiply all results by zero

         CALL SetErrStat(ErrID_Fatal, TRIM(p%OutParam(I)%Name)//" is not an available output channel.",ErrStat,ErrMsg,RoutineName)
      END IF

After this addition the outputs are updated properly, and I could see the FAST output for distributed forces in the local coordinates NDFl in OutData. This allows me to valid the transformation of the calculated Coriolis forces form local to global coordinates.

The figure above shows a good agreement between the calculated Coriolis forces in the local coordinates NCoriolisFland the transformed local distributed forces by FAST NDFl. Where, the transformed force in the global coordinates NCoriolisFg is the one added to u%DistrLoad%Force as below:

!! calculate the blade motion orientation tensor at each node 
   DO i=1,p%elem_total

      CALL BD_RotationalInterpQP( i, p, x, m )

      DO j=1,p%nqp

         CALL BD_CrvCompose(cc,m%qp%uuu(4:6,j,i),p%uu0(4:6,j,i),FLAG_R1R2)   ! cc = m%qp%uuu(4:6,j,i) composed with p%uu0(4:6,j,i)
         CALL BD_CrvCompose( cc0, p%Glb_crv, cc, FLAG_R1R2 )
         CALL BD_CrvMatrixR(cc,temp_R)   ! returns temp_R (the transpose of the DCM orientation matrix)

               ! Store the DCM for the j'th node of the i'th element (in FAST coordinate system)            
                BladeMotionOrientation(1:3,1:3,j) = TRANSPOSE(temp_R)
                BladeMotionOrientationTRANSPOSE(1:3,1:3,j) = TRANSPOSE(BladeMotionOrientation(1:3,1:3,j))               

               ! Applied distributed Coriolis forces at Node j  expressed in g, in N/m
                m%DistrCoriolisForceG( :,j ) = MATMUL(BladeMotionOrientationTRANSPOSE(:,:,j), m%DistrCoriolisForce( :,j))
                u%DistrLoad%Force(1:3,j)        =  m%DistrCoriolisForceG( 1:3,j ) 

      END DO
   END DO

However, in the example above the impact of the Coriolis force on the rotor speed can’t be seen. I would expect in the absence of aerodynamic and generator torque, the rotor would slow down between 10 and 50 s when a positive y component of Coriolis force is added, and the rotor would speed up between 60 and 100 s when a negative y component of Coriolis force is added. In this example I disabled the AeroDyn and ServoDyn and enable only the GenDOF.

Even when I apply a constant Coriolis force in the global coordinates, the rotor speed behavior doesn't seem as expected.

1. Do you see anything wrong in my implantation of the Coriolis force?

2. Do I have to change any other parameters in the simulation control ?

Best regards,

[jjonkman]

Dear @LaurenceWETI,

Your implementation now finally matches what I would expect to see. If I understand correctly, your results show that the force you've added is effecting the rotor speed, but it is not resulting in the rotor speed change that you expect. Are you applying your additional force in both BD_UpdateStates() and BD_CalcOutput()? Given the complicated nature of the force you've added, I can't say exactly how I'd expect the rotor speed to change, so, I would start by applying very simple force time histories to verify their effect.

Best regards,

[LaurenceWETI]

Dear @jjonkman,

Thank you very much for your quick reply.

That sounds good that the implementation now match what you would expect.

Regarding the simple force time histories would you please suggest a force time history that could slow down/ speed up the rotor speed? If I understand Figure 4.28, I would expect that positive force components in the yl direction should lead to slow down the rotor speed and vice versa. Therefore, I applied Coriolis forces in the previous example between 10 and 50 in the positive direction and between 60 and 100 in the negative direction. But I still can’t see the expected effect on the rotor speed.

Best regards,

[jjonkman]

Dear @LaurenceWETI,

I agree with your point, but considering that you are applying different loads to different nodes, this seems complicated to start with. I would first check the loads you are applying work well for the cantilevered beam case (no rotor rotation, by setting GenDOF = DrTrDOF = False, with RotSpeed = 0 rpm in ElastoDyn). And I would start with just one applied load at a given node. When doing this, does the blade deflection and reaction loads at the blade root make sense?

Best regards,

[LaurenceWETI]

Dear @jjonkman,

Thank you for your clarification.

As you suggested I applied just one force of 100 kN/m in yl at the blade tip and executed the test #26: NREL 5.0 MW Baseline Wind Turbine (Onshore); with the following settings in NRELOffshrBsline5MW_Onshore_ElastoDyn_BDoutputs.dat:

False         FlapDOF1 
False         FlapDOF2 
False         EdgeDOF  
False         TeetDOF  
False         DrTrDOF  
False         GenDOF   
False         YawDOF   
False         TwFADOF1 
False         TwFADOF2 
False         TwSSDOF1 
False         TwSSDOF2 
False         PtfmSgDOF
False         PtfmSwDOF
False         PtfmHvDOF
False         PtfmRDOF 
False         PtfmPDOF 
False         PtfmYDOF 

      0   RotSpeed

       0   PreCone(1)
       0   PreCone(2)
       0   PreCone(3)

        0   ShftTilt

The figure below shows the resulting tip translation deflections and root reaction force components.

1. Do the aforementioned results make sense for you?
2. Do you have any explanation why the results oscillate even when there is no rotation?

Best regards,

[jjonkman]

Dear @LaurenceWETI,

Overall, the results make sense to me. but the story is not fully clear because you are only showing results during the start-up transient. (The oscillations result because the blade is initialized with zero deflection, but the equilibrium condition with the applied load has deflection, resulting in a strong start-up transient. The oscillations should decay over time due to structural damping, but the mean value of the response is likely representative of the final static-equilibrium value.) The tip force in the y-direction results in deflection in the y-direction, which also leads to some deflection and reaction load in the axial direction.

Do the results make sense to you? A few clarifying questions:

• The tip load you are applying is expressed as a distributed force per unit length. So, the actual total force will depend on the element length. Can you confirm the element length and the total reaction load at the blade root.
• Does the natural frequency illustrated by the period of oscillation match what you expect?

Best regards,

[LaurenceWETI]

Dear @jjonkman,

Thank you for pointing that out.

After your explanation and trying to answer your questions the results make more since to me. Regarding your questions:

1. The length of the last element where the distributed force is applied is 0.15006 m. As the root reaction forces are expressed in r coordinates, which are illustrated in figure 2.29 in OpenFAST Documentation, I would expect a root reaction force of circa 15006 N in the opposite direction of the applied force. The figure below shows that the amplitude of the reaction force is as expected but the direction is not. Do I understand the r coordinates correctly?

2. Yes, the natural frequency illustrated by the period of oscillation (1.070 Hz) matches the 1st blade edgewise natural frequency of the NREL 5-MW wind turbine (1.0793 Hz).

• Is that possible to valid my implantation by proofing the law of conservation of momentum?

• Is the total momentum of the rotor calculated in the source code or implanted in the output data?

Best regards,

[jjonkman]

Dear @LaurenceWETI,

The signs of the root reaction loads are as I expect them. The reaction loads are computed as the load applied to the hub through the blade at the root. So, other than the small curvature and tip deflection, which changes the orientation of the "l" system relative to the "r" system, I would expect the signs of the forces in the "l" and "r" system to generally match.

Certainly conservation of momentum is preservered in the original BeamDyn code, but I don't believe that the total momentum is calculated directly within BeamDyn.

Best regards,