Implementation of external forces in BeamDyn

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

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,

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

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

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,

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

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 tower, 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

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,

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:

The TMD principle is based on a point mass that is attached to structure at a specific location. Whereas, the FW system deals with distributed mass along the structure. This means that the attached mass StC_M needs to be changed from a point mass to a distributed mass at specific nodes  along the structur.
The weight and the position of the distributed masses need to be dynamically and individually controlled for each blade within a given simulation.
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,

A few thoughts on this:

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?
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).
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

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

Dear @bjonkman, @jjonkman and @andrew-platt,
I did some changes in the source code of FAST_Library.f90, FAST_Registry.txt, FAST_Solver.f90 and Registry_BeamDyn.txt in order to pass 7 variables from Simulink through the S-Function to BeamDyn. These variables are still unused in BD. In the next figures I summarized the source code modifications:

FAST_Library.f90

FAST_Registry.txt

FAST_Solver.f90

Registry_BeamDyn.txt

The aforementioned files are attached to this message. I would appreciate it if someone could check the validity of the modified codes?
Modifications.zip
Matlab with the modified code could be compiled successfully, and the FAST_SFunc.mexw64 is generated. By running the open loop example in Matlab the simulation crushed with this error:

I couldn’t see anything strange in the .ech file that could refer to error.
Can anyone give me some hints to solve this error?
Best regards
Laurence

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 j and beam member i as 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 at j 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 FAT_Solver_ .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.

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

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

I would appreciate any help.
Best regards

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,

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

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,

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 FAT_Solver_ .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,

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,

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,

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,

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,

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 this 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,

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 INTENT on the model?
Best regards,

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,

Dear Jason,
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,

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,

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,

Dear @LaurenceWETI,
Can you share the full contents of SUBROUTINE BD_ChangeBladeMass()?
Best regards,
Best regards,

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,

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:

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?
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,

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,

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?

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,

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,

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,

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

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,

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,

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_ChangeBladeMass as 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,

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 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,

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 if condition 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  NCoriolisFl and 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.

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

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

Best regards,

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,

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,

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,

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.

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

Best regards,

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,

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:

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?

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,

Dear @jjonkman,
Thank for your help and your valuable tips.
In the same way I added and transformed local loads in BeamDyn coupled to FAST, I’m just trying to add local loads to BeamDyn stand-alone. Regarding the user documentation of OpenFAST the applied loads, including distributed, point, and tip-concentrated loads, are dead loads in the global coordinate system. Therefore, the additional local load must be converted to the global coordinate system before summing it with DistrLoad or PointLoad. The local coordinate system at each node depends on the initial position array, uu0, and the displacement (and rotation) array, uuu, at that node. These variables are outputs of BD_GaussPointData() subroutine. While these variables are not defined when I need to use them, I recalculated them in an additional subroutine BE_LocalLoads(), which is called in the very beginning of BD_CalcOutput() and BD_UpdateStates() subroutines.

!> @LaurenceWETI: This subroutine adds and transforms local loads.
SUBROUTINE BD_LocalLoads(p,x,m,u,ErrStat,ErrMsg)

   TYPE(BD_ParameterType),          INTENT(IN	)  :: p                 		!< Parameters
   TYPE(BD_ContinuousStateType), 	INTENT(IN   )  :: x           				!< Continuous states at t
   TYPE(BD_MiscVarType),            INTENT(INOUT)  :: m          				!< misc/optimization variables
   TYPE(BD_InputType),           	INTENT(INOUT)  :: u           				!< Inputs at t

   ! local variables
   REAL(BDKi)                  :: BladeMotionOrientation(3,3,p%ngp)
   REAL(BDKi)                  :: BladeMotionOrientationTRANSPOSE(3,3,p%ngp)
   REAL(BDKi)                  :: LocalForceL(3,p%node_total)
   REAL(BDKi)                  :: LocalForceG(3,p%node_total)

   REAL(BDKi)    			   :: hhx(p%node_elem)      						!< Shape function
   REAL(BDKi)    			   :: hpx(p%node_elem)      						!< Derivative of shape function
   REAL(BDKi)    			   :: uu0(6,p%ngp,p%elem_total)      				!< Initial position array at Gauss point
   REAL(BDKi)    			   :: E10(3)      									!< E10 = x_0^\prime at Gauss point
   REAL(BDKi)    			   :: uuu(p%dof_node,p%ngp,p%elem_total)      		!< Displacement(and rotation)  arrary at Gauss point
   REAL(BDKi)    			   :: uup(p%dof_node/2,p%ngp,p%elem_total)      	!< Derivative of displacement wrt axix at Gauss point
   REAL(BDKi)    			   :: E1(3)       									!< E1 = x_0^\prime + u^\prime at Gauss point
   INTEGER(IntKi)			   :: ErrStat     									!< Error status of the operation
   CHARACTER(*)  			   :: ErrMsg      									!< Error message if ErrStat /= ErrID_None

   REAL(BDKi)                  :: rrr(3)
   REAL(BDKi)                  :: rrp(3)
   REAL(BDKi)                  :: hhi
   REAL(BDKi)                  :: hpi
   REAL(BDKi)                  :: cc(3)
   REAL(BDKi)                  :: cc0(3)
   REAL(BDKi)                  :: rotu_temp(3)
   REAL(BDKi)                  :: rot_temp(3)
   REAL(BDKi)                  :: rot0_temp(3)
   REAL(BDKi)                  :: temp_R(3,3)
   INTEGER(IntKi)              :: inode
   INTEGER(IntKi)              :: temp_id
   INTEGER(IntKi)              :: temp_id2
   INTEGER(IntKi)              :: i
   INTEGER(IntKi)              :: j
   INTEGER(IntKi)              :: ErrStat2                     					! Temporary Error status
   INTEGER(IntKi)              :: nelem 										! number of elements
   INTEGER(IntKi)              :: igp											! number of Gauss point
   CHARACTER(ErrMsgLen)        :: ErrMsg2                      ! Temporary Error message
   CHARACTER(*), 					PARAMETER      :: RoutineName = 'BD_LocalLoads'

   ErrStat = ErrID_None
   ErrMsg  = ""

	BladeMotionOrientation			= 0.0_BDKi
	BladeMotionOrientationTRANSPOSE	= 0.0_BDKi
	uuu = 0.0_BDKi
	uup = 0.0_BDKi
	rrr = 0.0_BDKi
	rrp = 0.0_BDKi
	E10 = 0.0_BDKi
	uu0 = 0.0_BDKi

	DO nelem=1,p%elem_total

		DO igp=1,p%ngp

			DO inode=1,p%node_elem
				hhi = hhx(inode)
				hpi = hpx(inode)/p%Jacobian(igp,nelem)
				temp_id = (inode-1)*p%dof_node
				temp_id2 = (inode-1)*p%dof_node/2
				DO i=1,3
					uuu(i,igp,nelem) = uuu(i,igp,nelem) + hhi*m%Nuuu(temp_id+i)
					uup(i,igp,nelem) = uup(i,igp,nelem) + hpi*m%Nuuu(temp_id+i)
					rrr(i) = rrr(i) + hhi*m%Nrrr(temp_id2+i)
					rrp(i) = rrp(i) + hpi*m%Nrrr(temp_id2+i)
				ENDDO
			ENDDO

			E1 = 0.0_BDKi
			rotu_temp = 0.0_BDKi
			DO i=1,3
				E1(i) = E10(i) + uup(i,igp,nelem)
				rotu_temp(i) = m%Nuuu(i+3)
			ENDDO
			rot_temp(:)  = 0.0_BDKi
			rot0_temp(:) = 0.0_BDKi
			CALL BD_CrvCompose(rot_temp,rotu_temp,rrr,0,ErrStat2,ErrMsg2)
				CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
			DO i=1,3
				uuu(i+3,igp,nelem) = rot_temp(i)
				rot0_temp(i) = uu0(i+3,igp,nelem)
			ENDDO
		ENDDO
	ENDDO

               ! Define local point force at blade tip.
					LocalForceL(1,p%node_total) = 0.0D0
					LocalForceL(2,p%node_total) = 100000.0_BDKi
					LocalForceL(3,p%node_total) = 0.0D0

   DO i=1,p%elem_total

      DO j=1,p%ngp

         CALL BD_CrvCompose(cc,uuu(4:6,j,i),uu0(4:6,j,i),0,ErrStat2,ErrMsg2)   ! cc = uuu(4:6,j,i) composed with uu0(4:6,j,i)
			CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
         CALL BD_CrvCompose( cc0, p%Glb_crv, cc,0,ErrStat2,ErrMsg2)
			CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
         CALL BD_CrvMatrixR(cc,temp_R,ErrStat2,ErrMsg2)   ! 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
				LocalForceG( :,p%node_total) = MATMUL(BladeMotionOrientationTRANSPOSE(:,:,j), LocalForceL( :,p%node_total))
				u%PointLoad%Force(1:3,p%node_total)  =  LocalForceG( 1:3,p%node_total )

      END DO
   END DO

!			print*,"BD_LocalLoads"

END SUBROUTINE BD_LocalLoads

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

To verify the implementation is correct, I plotted the root reaction forces as response for a 100,000 N force applied between 100 and 200s at beam tip. The amplitude of the y-component of the root reaction force match the tip applied force.

Do my validation make sense?
Did I transform the forces correctly?
Why OutNd (Nodes whose values will be output) in BeamDyn stand-alone is limited to 6 outputs only?
How could I import external control signals to BeamDyn stand-alone in a similar way that I used through the S-Function in Simulink when BeamDyn coupled to FAST?

Best regards,

Dear @LaurenceWETI,
Regarding (1-2), I haven't reviewed your source code in detail (again, I'm not an expert on the BeamDyn source code), but what you describe makes sense to me.  I'm not sure I understand your verification test.  If you apply the load between 100-200s, what is happening at 50 s?
Regarding (3), there are a maximum of 9 output nodes (not 6) that can be set via the OutNd list and named outputs.  If you want more outputs than that (e.g., the full set of blade node motions and loads), use the nodal outputs instead: https://openfast.readthedocs.io/en/main/source/user/beamdyn/input_files.html#nodal-outputs.
Regarding (4), NREL has not developed a Simulink interface for any of the standalone module-level drivers, including BeamDyn.  I'm sure it is possible to create these, but you'd have to mimic features of the existing OpenFAST-Simulink interface and update the driver source code.
Best regards,

Dear @LaurenceWETI,
I'm not following.  I don't see the first block of code you show included in BeamDyn_IO.f90, where the calculation of nNodes is independent of the selected quadrature.  I so see the second block of code, where the calculation of nNodes depends on the selected quadrature.  Perhaps you are using an older version of the BeamDyn source code?
Best regards,

Dear @jjonkmann,
In the input file of BD stand-alone the root angular velocity, RootVel(4-6), can be set to constant values that do not change during the simulation. If the root motions need to be varied during the simulation, the source code must be modified. Do I understand that correctly?
Similarly to the initial rotor speed, RotSpeed, in ED input file, I need to run a simulation in BD stand-alone with an initial rotation speed, which can vary during the simulation. Is that possible in BD stand-alone? Can you give some hints who to do that?
I tried to couple my beam with OpenFAST, but when I run the simulation the solution does not converge after the maximum number of iterations. I’m not sure if this is because the beam properties that I’m using for the simulation in OpenFAST. Although, the simulation completed in BD stand-alone without any errors with the same beam properties. I can upload the input file if necessary?
Best regards,

Dear @LaurenceWETI,
Correct, the standalone BeamDyn driver is only set up to define a constant angular rotation of the root node about a hub axis.  If you want more complicated root motions, you must change the source code.  You would have to fully define the root motion (position, orientation, velocity, acceleration) as a function of time within the source code.  See the BeamDyn documentation for more information: https://openfast.readthedocs.io/en/main/source/user/beamdyn/input_files.html#blade-floating-reference-frame-parameters.
If your BeamDyn model solves as expected in the BeamDyn driver, but not in OpenFAST, this would suggest to me that you need to change the solver settings in OpenFAST, e.g., dropping the time step.  Due to the implicit loose coupling approach used by OpenFAST, models with ElastoDyn-BeamDyn coupling require smaller time steps than standalone BeamDyn models.
Best regards,

Dear @jjonkman,
Thank you for your quick response.
I would try first with BeamDyn coupled to OpenFAST bofore geting lost by the fully definition of the root motion :)
Yes, using smaller time steps did help. I set DT to 0.5 ms as suggested in this issue  #1300. But I still can't use the same beam properties that I used in the BeamDyn driver in OpenFAST. The beam that I solved in BeamDyn driver is very light and very stiff, e.g. m=0.015 kg/m and KShrFlp = 1.00E+17 Nm^2. Using these properties in OpenFAST leads again to convergence problem. Therefore, I had to change the masses and the stiffness to values that are similar to the values of 5MW reference rotor blade.  Do you think that such light and stiff beam can't be simulated in OpenFAST?
Best regards,

Dear @LaurenceWETI,
Those sectional properties do not sound physically correct.  And I would not recommend using BeamDyn for a very stiff blade.
OpenFAST should still be able to solve it, but would likely require a minuscule time step (too small to be practical).
Best regards,

Dear @jjonkman,
I set the following simulation parameters in the ElastoDyn_BDoutputs.dat in order to solve the aforementioned stiff BeamDyn model in OpenFAST:
False         FlapDOF1
False         FlapDOF2
False         EdgeDOF
False         TeetDOF
False         DrTrDOF
True          GenDOF
False         YawDOF
False         TwFADOF1
False         TwFADOF2
False         TwSSDOF1
False         TwSSDOF2
False         PtfmSgDOF
False         PtfmSwDOF
False         PtfmHvDOF
False         PtfmRDOF
False         PtfmPDOF
False         PtfmYDOF

0   OoPDefl
 0   IPDefl
 0   BlPitch(1)
 0   BlPitch(2)
 0   BlPitch(3)
 0   TeetDefl
90   Azimuth
 6   RotSpeed
 0   NacYaw
 0   TTDspFA
 0   TTDspSS
 0   PtfmSurge
 0   PtfmSway
 0   PtfmHeave
 0   PtfmRoll
 0   PtfmPitch
 0   PtfmYaw

       1   NumBl
      63   TipRad
     1.5   HubRad
       0	  PreCone(1)
       0	  PreCone(2)
       0	  PreCone(3)
       0   HubCM
       0   UndSling
       0   Delta3
       0   AzimB1Up
 -5.0191   OverHang
   1.912   ShftGagL
       0   ShftTilt
     1.9   NacCMxn
       0   NacCMyn
    1.75   NacCMzn
-3.09528   NcIMUxn
       0   NcIMUyn
 2.23336   NcIMUzn
 1.96256   Twr2Shft
    87.6   TowerHt
       0   TowerBsHt
       0   PtfmCMxt
       0   PtfmCMyt
       0   PtfmCMzt
       0   PtfmRefzt

OpenFAST requires a time step of 0.2 ms to solve my beam model.
I’m just wondering, the 1p excitation in the rotor speed signal where does it come from?

Best regards,

Dear everyone,
Dear @jjonkman,
is it possible to attaches the beam to the hub as pinned support, without the effect of the attached inertia of the drivetrain and the generator?
I tried to enable only the GenDOF, but I can't set the GenIner to zero (see above).
the problem is with 0.3 generator inertia the initial rotor speed didn't start from initial value of 10 rpm, it has an offset of 0.2 rpm. could anyone help me to simulate a free spinning beam in OpenFAST/Simulink.
Best regards,

Dear @LaurenceWETI,
In your OpenFAST model, the rigid-body rotational DOF of the drivetrain is included in ElastoDyn, and because the blades are not included in the ElastoDyn model when BeamDyn is enabled, ElastoDyn needs some nonzero-inertia associated with the with this rotational DOF.  You'll need to specify some physically realistic drivetrain rotational inertia via HubIner and/or GenIner.
Best regards,

Dear @jjonkman,
Thank you very much for your reply.
It looks now much better as I used the drivetrain rotational inertia of the NREL 5MW. But I still have a problem with the support between the beam and the hub. I would simulate a pinned support between the beam and the hub. In the figure bellows I plotted the reaction forces and moments at the blade root. I would expect that the root reaction moment should be almost by zero for a pinned support, which can’t be seen in my results.

Is it possible to model such kind of supports in OpenFAST/BeamDyn?

Best regards,

Dear @LaurenceWETI,
The ElastoDyn-BeamDyn coupling available in OpenFAST only supports a cantilevered connection between the blade (in BeamDyn) and hub (in ElastoDyn).
Best regards,

Dear everyone,
Dear @jjonkman,
I have a question regarding the calculation of the element length in BeamDyn.
My curved bema has the following geometry:

---------------------- GEOMETRY PARAMETER --------------------------------------
          1   member_total    - Total number of members (-)
         15   kp_total        - Total number of key points (-) [must be at least 3]
     1     15                 - Member number; Number of key points in this member
   kp_xr         kp_yr         kp_zr        initial_twist
   (m)            (m)          (m)            (deg)
0.000000E+00	0.000000E+00	0.000000E+00	0.000000E+00
1.000000E-02	0.000000E+00	1.000000E+00	0.000000E+00
4.000000E-02	0.000000E+00	2.000000E+00	0.000000E+00
9.000000E-02	0.000000E+00	3.000000E+00	0.000000E+00
1.600000E-01	0.000000E+00	4.000000E+00	0.000000E+00
2.498000E-01	0.000000E+00	4.999000E+00	0.000000E+00
2.500000E-01	0.000000E+00	5.000100E+00	0.000000E+00
5.625000E-01	0.000000E+00	7.500000E+00	0.000000E+00
9.980000E-01	0.000000E+00	9.999000E+00	0.000000E+00
1.000000E+00	0.000000E+00	1.000100E+01	0.000000E+00
1.210000E+00	0.000000E+00	1.100000E+01	0.000000E+00
1.440000E+00	0.000000E+00	1.200000E+01	0.000000E+00
1.690000E+00	0.000000E+00	1.300000E+01	0.000000E+00
1.960000E+00	0.000000E+00	1.400000E+01	0.000000E+00
2.250000E+00	0.000000E+00	1.500000E+01	0.000000E+00

Regarding my hand calculation the distance along the local axis between keypoints (kp) 1-5, 6-9 and 10-15 should be 5.007241985, 5.05747081 and 5.156729137, respectively. However, BeamDyn results different distances between the aforementioned kps.
Root_totoal_length =
 4.9435
Body_totoal_length =
 5.0415
Tip_totoal_length =
 5.2384

The code that I used to calculate and display the distances in MATLAB is as follow:
DO i=1,p%elem_total !
      DO j=1,p%nqp!

   		IF (j >= 1.AND. j <= 5) THEN
   			Root_length(i) = Root_length(i) + p%QPtWeight(j)*p%Jacobian(j,i)
   		ELSEIF (j > 5 .AND. j < 10) THEN
   			Body_length(i) = Body_length(i) + p%QPtWeight(j)*p%Jacobian(j,i)
   		ELSEIF (j >= 10 .AND. j <= 15) THEN
   			Tip_length(i) = Tip_length(i) + p%QPtWeight(j)*p%Jacobian(j,i)
   		ELSE
   		ENDIF
      ENDDO
      p%Root_totoal_length	= p%Root_totoal_length	+ Root_length(i)
      p%Body_totoal_length 	= p%Body_totoal_length 	+ Body_length(i)
      p%Tip_totoal_length	= p%Tip_totoal_length 	+ Tip_length(i)
  ENDDO
   	CALL WrScr("Root_totoal_length =")
   	CALL WrScr( trim( Num2Lstr( p%Root_totoal_length ) ) )
   	CALL WrScr("Body_totoal_length =")
   	CALL WrScr( trim( Num2Lstr( p%Body_totoal_length ) ) )
   	CALL WrScr("Tip_totoal_length =")
   	CALL WrScr( trim( Num2Lstr( p%Tip_totoal_length ) ) )

Do you have any idea why the distances in BeamDyn differ from my hand calculation ?

Best regards

Dear @LaurenceWETI,
The keypoint line does not explicitly define the location of the quadrature points.  Rather the keypoint line defines the line that the quadrature points fall on.  The actual locations of the quadrature points is from the Eta values for each blade input station in the blade definition file (the single values preceding each stiffness and mass matrix set).  These Eta values are placed along the curve of the keypoint line so the distance between quadrature points is along the curve.
Andy

Dear @andrew-platt,
Thank you for your quick reply.
I derived the local beam reference axis and the nondimensional parameter Eta, η, in the beam definition file based on the keypoint line as follow:
Hand calculation:
kp | distance along local axis | Eta |   |   |   |
-- | -- | -- | -- | -- | -- | --
(-) | (m) | η(-) |   |   |   |
1 | 0.00E+00 | 0.000E+00 |   |   |   |
2 | 1.00E+00 | 6.570E-02 |   |   |   |
3 | 2.00E+00 | 1.314E-01 |   |   |   |
4 | 3.00E+00 | 1.972E-01 |   |   |   |
5 | 4.00E+00 | 2.631E-01 |   |   |   |
6 | 5.01E+00 | 3.290E-01 | kp1-6 | Root_total_lengeh= | 5.01E+00
7 | 5.01E+00 | 3.291E-01 |   |   |   |
8 | 7.53E+00 | 4.945E-01 |   |   |   |
9 | 1.01E+01 | 6.612E-01 | kp9-6 | Body_total_lengeh= | 5.06E+00
10 | 1.01E+01 | 6.614E-01 |   |   |   |
11 | 1.11E+01 | 7.284E-01 |   |   |   |
12 | 1.21E+01 | 7.958E-01 |   |   |   |
13 | 1.31E+01 | 8.635E-01 |   |   |   |
14 | 1.42E+01 | 9.316E-01 |   |   |   |
15 | 1.52E+01 | 1.000E+00 | kp15-9 | Tip_total_lengeh= | 5.16E+00
BeamDyn calculation:
 DO i=1,p%elem_total !
       DO j=1,p%nqp!

			IF (j >= 1.AND. j <= 6) THEN
				Root_length(i) = Root_length(i) + p%QPtWeight(j)*p%Jacobian(j,i)
			ELSEIF (j > 6 .AND. j < 10) THEN
				Body_length(i) = Body_length(i) + p%QPtWeight(j)*p%Jacobian(j,i)
			ELSEIF (j >= 10 .AND. j <= 15) THEN
				Tip_length(i) = Tip_length(i) + p%QPtWeight(j)*p%Jacobian(j,i)
			ELSE
			ENDIF
       ENDDO
       p%Root_totoal_length	= p%Root_totoal_length	+ Root_length(i)
       p%Body_totoal_length 	= p%Body_totoal_length 		+ Body_length(i)
       p%Tip_totoal_length	= p%Tip_totoal_length 	+ Tip_length(i)
   ENDDO
		CALL WrScr("Root_totoal_length =")
		CALL WrScr( trim( Num2Lstr( p%Root_totoal_length ) ) )
		CALL WrScr("Body_totoal_length =")
		CALL WrScr( trim( Num2Lstr( p%Body_totoal_length ) ) )
		CALL WrScr("Tip_totoal_length =")
		CALL WrScr( trim( Num2Lstr( p%Tip_totoal_length ) ) )

Root_totoal_length =
 4.9435
Body_totoal_length =
 5.0415
Tip_totoal_length =
 5.2384

Do I use the correct parameters,  p%QPtWeight(j) and  p%Jacobian(j,i), to calculate the element length ?

Is there any other method to derive the Eta values?

Can you see any mistakes in my calculation?

Best regards,

Dear everyone,
Dear @jjonkman and @andrew-platt
I added external force at several nodes of the blades as distributed forces. In the figure below you can see the additional loads and the corresponding applied distributed forces expressed in the local coordinates.

As you see, the applied distributed forces B1N4DFyl for the first blade are matching the additional external distributed forces B1N4AddFyl. However, they deviate for second and the third blades.
Here is how I transformed the forces from local to global coordinates.
!!calculate the blade motion orientation tensor at each node, and transform the local external forces from local to global coordinates.
   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 external forces at Node j  expressed in g, in N/m
				m%DistrExternalForceG( :,j )	= MATMUL(BladeMotionOrientationTRANSPOSE(:,:,j), m%DistrExternalForceL(:,j))

               ! Applied distributed external forces at Node j  expressed in g, in N/m
				u%DistrLoad%Force(1:3,j)		=  m%DistrExternalForceG( 1:3,j )
      END DO
   END DO

Do you have any idea why the NDFl are deviating from the NADDFL for blades 2 and 3?
Best regards,

I can't spot the error.

@LaurenceWETI, can you plot the values in m%DistrExternalForceL?  I also don't see any issues with the code above.

Dear @jjonkman and @andrew-platt
thank you very much for your replies.
The figures below compare the x, y, and z components of the external applied distributed forces, m%DistrExternalForceL, and the applied distributed forces, NDFl,expressed in the local coordinates for each of the blades.

I think the differences in the y and z components of blades 2 and 3 should come from the transformation of, u%DistrLoad%Force,from global to local coordinates. Do you agree?

I can’t also understand why do the x components of, NDFxl, have the sine wave form for all blades?

Best regards

Dear @LaurenceWETI,
Where are you adding the calculations you shared in your post from 3 days ago within the BeamDyn source code?  Are you running BeamDyn within OpenFAST?  Are there other loads being applied to the blades?  Is the rotor spinning?
Best regards,

Dear @jjonkman,
Thank you for your reply.

The calculation is added at the very beginning of BD_UpdateStates and BD_CalcOutput subroutines

Yes, BeamDyn is running within OpenFAST Test # 26 with deactivated AeroDyn and ServoDyn . The Gravity equals zero.

All external loads are summed up in m%DistrExternalForceL. There are no other external loads.

Yes, the rotor is spinning. The rotational speed equals 12.1 rpm

Best regards

Dear @LaurenceWETI,
I'm not sure what is going on.  I would suggest stepping through with the debugger.  The sine wave seems to have the period equivalent of the rotor speed, but I'm not sure why the output in local coordinates is not the same as you are specifying.
Best regards,

Dear @LaurenceWETI,
The values for p%uu0() are parameters that should remain constant thoughout the entire BeamDyn simulation.  As you noted, they are identical for all blades.
During initialization, the root reference orientation is passed into BeamDyn (this is the root orientation of a given with no initial azimuth rotation for the rotor).  The blade reference frame defined by p%uu0 is defined relative to the initial root orientation for each blade.  Thus it doesn't matter that blades 2 and 3 are pointed in different directions since p%uu0 is relative to their respective root reference frames, and therefore p%uu0 will be identical for all blades.
It is possible that this is causing issues with the orientations of your applied loads on blades 2 and 3.  I don't know exactly what you've changed in the code to apply the loads, but I am suspicious that you need to include the root reference orientation in your load application.  This is stored in the variable p%GlbRot and gives the orientation of the root reference frame.  You might be missing a transformation like m%DistrExternalForceL = matmul(transpose(p%GlbRot,<ForceInGlobalCoords>)) to transform from global to the blade reference coordinates, then apply whatever other transforms you need (I'm not sure which coordinate frame you are applying your loads in -- also not sure which order of orientation changes you will need).
On the plots in the above comment, I am unsure what the x-axis is.  Do I assume correctly that it is time?
Hopefully this is helpful in finding the issue. Regards,
Andy

Dear @andrew-platt,
Thank you very much for pointing that out.
I thought that p%uu0 defines the blade reference frame. That is why I expected to see the azimuth in its value for blade 2 and 3. Where I can see the definition of the blade reference frame in the source code of BeamDyn?
May I briefly explain what I’ve changed in the code to apply the loads:
The calculation of the distributed external forces m%DistrExternalFoeceL is based on the local components of the rotational velocity at the quadrature points  m%qp%vvv(4:6,:,:)  :
m%DistrExternalForceL(:,:) = u%ExternalVariable * m%qp%vvv(4:6,:,:)
wehre u%ExternalVariable is an external variable that is imported through the FAST S-Func and called for each blade at the end of BD_InputSolve subroutine in FAST_Solver.f90 as follow:
DO K = 1,p_FAST%nBeams ! Loop through all blades

	IF ( K == 1) THEN
				BD%Input(1,k)% ExternalVariable       	=  m_FAST%ExternInput% ExternalVariableB1
	ELSEIF ( K == 2) THEN
				BD%Input(1,k)% ExternalVariable       	=  m_FAST%ExternInput% ExternalVariableB2
	ELSE
				BD%Input(1,k)% ExternalVariable       	=  m_FAST%ExternInput% ExternalVariableB3
END DO

The calculation of the external forces takes place in the additional subroutine BD_ExternalForces. This subroutine is called at the very beginning of BD_UpdateStates and BD_CalcOutput subroutines. Since, the distributed applied forces u%DistrLoad%Force are expressed in the global coordinates when BD_ExternalForces subroutines is called, the external forces must be transformed form the local to the global coordinates. The code used for the transformation from local to global coordinates is posted in the my post form the last week:
You mentioned that I need to use p%GlbRot to transform my external forces from global to local coordinates. Do you mean that I need to calculate my external force in the global coordinates and then add them to u%DistrLoad%Force. Are  y%BldMotion%RotationVel the transformed rotational velocity in FAST coordinate system?  if I replace y%BldMotion%RotationVel with m%qp%vvv in the equation above, the external forces would be then expressed in the global coordinates? and can be added to u%DistrLoad%Force without any transformation? I would try this and let you know.
But I still wondering why the transformation from local to global coordinates works for blade 1 and doesn’t for blades 2 and 3? Do you see anything wrong in my implementation?
Best regards,

Dear @andrew-platt, Dear @jjonkman,
In the example below I just tried to simplify the calculation of the local external loads to a given constant distributed force of 20 N/m in a given period of 5 s between the 5th and 10th seconds for a total simulation time of 15 s. The local loads are transformed in an additional subroutine named BD_TransfLocal2Global, which is called in the beginning of BD_UpdateStates and BD_CalcOutput subroutines. To insure that this is the only modification that I’ve done in the source code of BeamDyn.f90, I cloned a new version of OpenFAST and implemented this simple example inside.
These are the sections where the code is modified:
SUBROUTINE BD_UpdateStates( t, n, u, utimes, p, x, xd, z, OtherState, 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
   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;
                                                                !!   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_TransfLocal2Global.
   INTEGER(IntKi)                               :: i           ! generic loop counter
   DO i=1,size(utimes)
		IF( t > 5 .AND. t < 10 ) THEN
			CALL BD_TransfLocal2Global(p,x,m,u(i),ErrStat,ErrMsg)
		ELSE
			u(i)%DistrLoad%Force(:,:)		= 0.0_ReKi
		ENDIF
   ENDDO

   ! Initialize ErrStat

SUBROUTINE BD_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, NeedWriteOutput )

   REAL(DbKi),                   INTENT(IN   )  :: t           !< Current simulation time in seconds
   TYPE(BD_InputType),           INTENT(INOUT)  :: u           !< Inputs at t
   TYPE(BD_ParameterType),       INTENT(IN   )  :: p           !< Parameters
   TYPE(BD_ContinuousStateType), INTENT(IN   )  :: x           !< Continuous states at t
   TYPE(BD_DiscreteStateType),   INTENT(IN   )  :: xd          !< Discrete states at t
   TYPE(BD_ConstraintStateType), INTENT(IN   )  :: z           !< Constraint states at t
   TYPE(BD_OtherStateType),      INTENT(IN   )  :: OtherState  !< Other states at t
   TYPE(BD_OutputType),          INTENT(INOUT)  :: y           !< Outputs computed at t (Input only so that mesh con-
                                                               !!   nectivity information does not have to be recalculated)
   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
   LOGICAL,          OPTIONAL,   INTENT(IN   )  :: NeedWriteOutput     !< Flag to determine if WriteOutput values need to be calculated in this call

   TYPE(BD_ContinuousStateType)                 :: x_tmp
   TYPE(BD_OtherStateType)                      :: OtherState_tmp
   INTEGER(IntKi)                               :: i           ! generic loop counter
   INTEGER(IntKi)                               :: nelem       ! loop over elements
   REAL(ReKi)                                   :: AllOuts(0:MaxOutPts)
   INTEGER(IntKi)                               :: ErrStat2                     ! Temporary Error status
   CHARACTER(ErrMsgLen)                         :: ErrMsg2                      ! Temporary Error message
   CHARACTER(*), PARAMETER                      :: RoutineName = 'BD_CalcOutput'
   LOGICAL                                      :: CalcWriteOutput

!@LaurenceWETI: Call Subroutine BD_TransfLocal2Global.
		IF( t > 5 .AND. t < 10 ) THEN
      CALL BD_TransfLocal2Global(p,x,m,u,ErrStat,ErrMsg)
  ELSE
	u%DistrLoad%Force				= 0.0_ReKi
   ENDIF

   ! Initialize ErrStat

END SUBROUTINE BD_ComputeElementMass

!-----------------------------------------------------------------------------------------------------------------------------------
!> @LaurenceWETI: This subroutine transforms loads in local coordinates to global coordinates.
SUBROUTINE BD_TransfLocal2Global(p,x,m,u,ErrStat,ErrMsg)

   TYPE(BD_ParameterType),          INTENT(IN   )  :: p                 !< Parameters
   TYPE(BD_ContinuousStateType), 	INTENT(IN   )  :: x           		!< Continuous states at t
   TYPE(BD_MiscVarType),            INTENT(INOUT)  :: m          		!< misc/optimization variables
   TYPE(BD_InputType),           	INTENT(INOUT)  :: u           		!< Inputs at t
																		!!   nectivity information does not have to be recalculated)

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

   CHARACTER(*), 					PARAMETER      :: RoutineName = 'BD_TransfLocal2Global'

   ! local variables
   INTEGER(IntKi)              :: i
   INTEGER(IntKi)              :: j
   REAL(BDKi)                  :: cc(3)
   REAL(BDKi)                  :: cc0(3)
   REAL(BDKi)                  :: temp_R(3,3)
   REAL(BDKi)                  :: BladeMotionOrientation(3,3,p%nqp)
   REAL(BDKi)                  :: BladeMotionOrientationTRANSPOSE(3,3,p%nqp)
   REAL(BDKi)                  :: DistrExternalForceL(3,p%nqp)				!< fluid momentum forces in local coordinates at Element start
   REAL(BDKi)                  :: DistrExternalForceG(3,p%nqp)				!< fluid momentum forces in local coordinates at Element end

   ErrStat = ErrID_None
   ErrMsg  = ""

	DistrExternalForceL				= 0.0_ReKi
	DistrExternalForceG				= 0.0_ReKi

!_____________________________________________________________________________
!!calculate the blade motion orientation tensor at each node, and transform the local forces to global coordinates.
   DO i=1,p%elem_total

      ! Set the initial displacements: p%uu0, p%rrN0, p%E10
      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))

				DistrExternalForceL(:,j) = 20.00
               ! Applied distributed Coriolis and inertial forces at Node j  expressed in g, in N/m
				DistrExternalForceG( :,j )	= MATMUL(BladeMotionOrientationTRANSPOSE(:,:,j), DistrExternalForceL(:,j))

               ! Applied distributed external forces at Node j  expressed in g, in N/m
				u%DistrLoad%Force(:,j)		= DistrExternalForceG(:,j )
      END DO
   END DO

END SUBROUTINE BD_TransfLocal2Global

!-----------------------------------------------------------------------------------------------------------------------------------
!> This routine alters the RootMotion inputs based on the pitch-actuator parameters and discrete states
SUBROUTINE PitchActuator_SetBC(p, u, xd, AllOuts)

The figure below shows the applied distributed loads at nod 4 (NREL5MW blade station 20) expressed in the local coordinates for blade 1, 2, and 3.  The simulation scenario used here is the 5MW_Land_BD_DLL_WTurb with deactivated gravity, servodyn, aerody, infolow, and all degrees of freedom in elastodyn. The rotor is not spinning (RotSpeed=0).

It is obvious that the applied distributed loads for blades 2 and 3 are not matching with the given constant distributed force of 20 N/m.
With this example the bug can be excluded in the code that I used for the transformation from the local to the global coordinates.
Do you agree with my conclusion? Or do see anything unclear in my example?
Best regards,
PS: I attached also the BeamDyn.f90 file that I modified.
BeamDyn.txt


[rocfft] Include  as std::cerr is referenced.

resolve TC build
the issue is caused by missing header, 'std::cerr' is referenced but  is not included

Hi mhbliao,
Any PR committed to develop branch has to be:

pass CI tests
get at least one peer review approve.

Please follow the rules and try it again if you really have something to commit.

http://ontrack-internal.amd.com/browse/SWDEV-193409

I have submitted another pull request. Please review that simple test as soon as possible as rocFFT failed on TC builds for a while.
Thanks- Michael
On Thursday, June 20, 2019, 2:24:33 PM EDT, feizheng10<EMAIL_ADDRESS>wrote:

Hi mhbliao,
Any PR committed to develop branch has to be:

pass CI tests
get at least one peer review approve.

Please follow the rules and try it again if you really have something to commit.
—
You are receiving this because you modified the open/close state.
Reply to this email directly, view it on GitHub, or mute the thread.

BTW, how often you merge develop into master, we also need to merge that simple change into master ASAP to unlock TC builds.
Thanks- Michael
On Thursday, June 20, 2019, 2:59:45 PM EDT, Michael Liao<EMAIL_ADDRESS>wrote:

I have submitted another pull request. Please review that simple test as soon as possible as rocFFT failed on TC builds for a while.
Thanks- Michael
On Thursday, June 20, 2019, 2:24:33 PM EDT, feizheng10<EMAIL_ADDRESS>wrote:

Hi mhbliao,
Any PR committed to develop branch has to be:

pass CI tests
get at least one peer review approve.

Please follow the rules and try it again if you really have something to commit.
—
You are receiving this because you modified the open/close state.
Reply to this email directly, view it on GitHub, or mute the thread.