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): #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 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.

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

[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 INTENT on 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,

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

[LaurenceWETI]

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.

1. Do my validation make sense?
2. Did I transform the forces correctly?
3. Why OutNd (Nodes whose values will be output) in BeamDyn stand-alone is limited to 6 outputs only?
4. 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,

[jjonkman]

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,

[LaurenceWETI]

Dear @jjonkman,

Thank you very much for your quick reply.
Regarding 1: sorry it was my fault, the local force at blade tip is applied starting from 50 s instead of 100 s.
Regarding 3: this is what I would expect, but when I set OutNd to values above 6, I’ve got this error massage:

This error message comes from this check in BeamDyn_IO.f90:

  ! Check to see if all OutNd(:) analysis points are existing analysis points:
      nNodes = (InputFileData%order_elem + 1)*InputFileData%member_total  ! = p%node_elem*p%elem_total (number of nodes on y%BldMotion mesh)
       
      do j=1,InputFileData%NNodeOuts
         if ( InputFileData%OutNd(j) < 1_IntKi .OR. InputFileData%OutNd(j) > nNodes ) then
            call SetErrStat( ErrID_Fatal, ' All OutNd values must be between 1 and '//&
                    trim( Num2LStr( nNodes ) )//' (inclusive).', ErrStat, ErrMsg, RoutineName )
            exit ! stop checking this loop
         end if
      end do

Where order_elem is set to 5 in the input file.
This check in BeamDyn coupled to FAST is extended to accept the values of blade stations as follow:

   ! Check to see if all OutNd(:) analysis points are existing analysis points:
!bjj: FIXME!!!! THIS ISN'T NECESSARIALLY THE NUMBER OF NODES ON THIS MESH
!      IF (p%BldMotionNodeLoc == BD_MESH_QP) THEN
         if(InputFileData%quadrature .EQ. TRAP_QUADRATURE) then
            nNodes = (InputFileData%InpBl%station_total - 1)*InputFileData%refine + 1  ! number of nodes on y%BldMotion mesh
         else
            nNodes = (InputFileData%order_elem + 1)*InputFileData%member_total  ! = p%nodes_per_elem*p%elem_total (number of nodes on y%BldMotion mesh)
         end if

Do I have to change this check also as described in BeamDyn coupled to FAST?

Best regards,

[jjonkman]

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,

[LaurenceWETI]

Dear @jjonkman,

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,

[jjonkman]

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,

[LaurenceWETI]

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,

[jjonkman]

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,

[LaurenceWETI]

Dear @jjonkman,

Since, the BeamDyn driver is only set up to define a constant angular rotation of the root node about a hub axis, I’m trying to solve my stiff beam in OpenFAST/BeamDyn. I set the following simulation parameters in ElastoDyn_BDoutputs.dat in order to solve my model for free rotational speed.

------- ELASTODYN v1.03.* INPUT FILE -------------------------------------------
NREL 5.0 MW Baseline Wind Turbine for Use in Offshore Analysis. Properties from Dutch Offshore Wind Energy Converter (DOWEC) 6MW Pre-Design (10046_009.pdf) and REpower 5M 5MW (5m_uk.pdf)
---------------------- SIMULATION CONTROL --------------------------------------
False         Echo        - Echo input data to "<RootName>.ech" (flag)
          3   Method      - Integration method: {1: RK4, 2: AB4, or 3: ABM4} (-)
"DEFAULT"     DT          - Integration time step (s)
---------------------- DEGREES OF FREEDOM --------------------------------------
False         FlapDOF1    - First flapwise blade mode DOF (flag)
False         FlapDOF2    - Second flapwise blade mode DOF (flag)
False         EdgeDOF     - First edgewise blade mode DOF (flag)
False         TeetDOF     - Rotor-teeter DOF (flag) [unused for 3 blades]
False         DrTrDOF     - Drivetrain rotational-flexibility DOF (flag)
True          GenDOF      - Generator DOF (flag)
False         YawDOF      - Yaw DOF (flag)
False         TwFADOF1    - First fore-aft tower bending-mode DOF (flag)
False         TwFADOF2    - Second fore-aft tower bending-mode DOF (flag)
False         TwSSDOF1    - First side-to-side tower bending-mode DOF (flag)
False         TwSSDOF2    - Second side-to-side tower bending-mode DOF (flag)
False         PtfmSgDOF   - Platform horizontal surge translation DOF (flag)
False         PtfmSwDOF   - Platform horizontal sway translation DOF (flag)
False         PtfmHvDOF   - Platform vertical heave translation DOF (flag)
False         PtfmRDOF    - Platform roll tilt rotation DOF (flag)
False         PtfmPDOF    - Platform pitch tilt rotation DOF (flag)
False         PtfmYDOF    - Platform yaw rotation DOF (flag)
---------------------- INITIAL CONDITIONS --------------------------------------
          0   OoPDefl     - Initial out-of-plane blade-tip displacement (meters)
          0   IPDefl      - Initial in-plane blade-tip deflection (meters)
          0   BlPitch(1)  - Blade 1 initial pitch (degrees)
          0   BlPitch(2)  - Blade 2 initial pitch (degrees)
          0   BlPitch(3)  - Blade 3 initial pitch (degrees) [unused for 2 blades]
          0   TeetDefl    - Initial or fixed teeter angle (degrees) [unused for 3 blades]
         90   Azimuth     - Initial azimuth angle for blade 1 (degrees)
		 10   RotSpeed    - Initial or fixed rotor speed (rpm)
          0   NacYaw      - Initial or fixed nacelle-yaw angle (degrees)
          0   TTDspFA     - Initial fore-aft tower-top displacement (meters)
          0   TTDspSS     - Initial side-to-side tower-top displacement (meters)
          0   PtfmSurge   - Initial or fixed horizontal surge translational displacement of platform (meters)
          0   PtfmSway    - Initial or fixed horizontal sway translational displacement of platform (meters)
          0   PtfmHeave   - Initial or fixed vertical heave translational displacement of platform (meters)
          0   PtfmRoll    - Initial or fixed roll tilt rotational displacement of platform (degrees)
          0   PtfmPitch   - Initial or fixed pitch tilt rotational displacement of platform (degrees)
          0   PtfmYaw     - Initial or fixed yaw rotational displacement of platform (degrees)
---------------------- TURBINE CONFIGURATION -----------------------------------
          1   NumBl       - Number of blades (-)
         15   TipRad      - The distance from the rotor apex to the blade tip (meters)
		  0   HubRad      - The distance from the rotor apex to the blade root (meters)
          0	  PreCone(1)  - Blade 1 cone angle (degrees)
          0	  PreCone(2)  - Blade 2 cone angle (degrees)
          0	  PreCone(3)  - Blade 3 cone angle (degrees) [unused for 2 blades]
          0   HubCM       - Distance from rotor apex to hub mass [positive downwind] (meters)
          0   UndSling    - Undersling length [distance from teeter pin to the rotor apex] (meters) [unused for 3 blades]
          0   Delta3      - Delta-3 angle for teetering rotors (degrees) [unused for 3 blades]
          0   AzimB1Up    - Azimuth value to use for I/O when blade 1 points up (degrees)
    -5.0191   OverHang    - Distance from yaw axis to rotor apex [3 blades] or teeter pin [2 blades] (meters)
      1.912   ShftGagL    - Distance from rotor apex [3 blades] or teeter pin [2 blades] to shaft strain gages [positive for upwind rotors] (meters)
          0   ShftTilt    - Rotor shaft tilt angle (degrees)
        1.9   NacCMxn     - Downwind distance from the tower-top to the nacelle CM (meters)
          0   NacCMyn     - Lateral  distance from the tower-top to the nacelle CM (meters)
       1.75   NacCMzn     - Vertical distance from the tower-top to the nacelle CM (meters)
   -3.09528   NcIMUxn     - Downwind distance from the tower-top to the nacelle IMU (meters)
          0   NcIMUyn     - Lateral  distance from the tower-top to the nacelle IMU (meters)
    2.23336   NcIMUzn     - Vertical distance from the tower-top to the nacelle IMU (meters)
    1.96256   Twr2Shft    - Vertical distance from the tower-top to the rotor shaft (meters)
       87.6   TowerHt     - Height of tower above ground level [onshore] or MSL [offshore] (meters)
          0   TowerBsHt   - Height of tower base above ground level [onshore] or MSL [offshore] (meters)
          0   PtfmCMxt    - Downwind distance from the ground level [onshore] or MSL [offshore] to the platform CM (meters)
          0   PtfmCMyt    - Lateral distance from the ground level [onshore] or MSL [offshore] to the platform CM (meters)
          0   PtfmCMzt    - Vertical distance from the ground level [onshore] or MSL [offshore] to the platform CM (meters)
          0   PtfmRefzt   - Vertical distance from the ground level [onshore] or MSL [offshore] to the platform reference point (meters)
---------------------- MASS AND INERTIA ----------------------------------------
          0   TipMass(1)  - Tip-brake mass, blade 1 (kg)
          0   TipMass(2)  - Tip-brake mass, blade 2 (kg)
          0   TipMass(3)  - Tip-brake mass, blade 3 (kg) [unused for 2 blades]
		  0   HubMass     - Hub mass (kg)
		  0   HubIner     - Hub inertia about rotor axis [3 blades] or teeter axis [2 blades] (kg m^2)
    0.3   GenIner     - Generator inertia about HSS (kg m^2)
		  0   NacMass     - Nacelle mass (kg)
		  0   NacYIner    - Nacelle inertia about yaw axis (kg m^2)
          0   YawBrMass   - Yaw bearing mass (kg)
          0   PtfmMass    - Platform mass (kg)
          0   PtfmRIner   - Platform inertia for roll tilt rotation about the platform CM (kg m^2)
          0   PtfmPIner   - Platform inertia for pitch tilt rotation about the platform CM (kg m^2)
          0   PtfmYIner   - Platform inertia for yaw rotation about the platform CM (kg m^2)
---------------------- BLADE ---------------------------------------------------
         17   BldNodes    - Number of blade nodes (per blade) used for analysis (-)
"NRELOffshrBsline5MW_Blade.dat"    BldFile(1)  - Name of file containing properties for blade 1 (quoted string)
"NRELOffshrBsline5MW_Blade.dat"    BldFile(2)  - Name of file containing properties for blade 2 (quoted string)
"NRELOffshrBsline5MW_Blade.dat"    BldFile(3)  - Name of file containing properties for blade 3 (quoted string) [unused for 2 blades]
---------------------- ROTOR-TEETER --------------------------------------------
          0   TeetMod     - Rotor-teeter spring/damper model {0: none, 1: standard, 2: user-defined from routine UserTeet} (switch) [unused for 3 blades]
          0   TeetDmpP    - Rotor-teeter damper position (degrees) [used only for 2 blades and when TeetMod=1]
          0   TeetDmp     - Rotor-teeter damping constant (N-m/(rad/s)) [used only for 2 blades and when TeetMod=1]
          0   TeetCDmp    - Rotor-teeter rate-independent Coulomb-damping moment (N-m) [used only for 2 blades and when TeetMod=1]
          0   TeetSStP    - Rotor-teeter soft-stop position (degrees) [used only for 2 blades and when TeetMod=1]
          0   TeetHStP    - Rotor-teeter hard-stop position (degrees) [used only for 2 blades and when TeetMod=1]
          0   TeetSSSp    - Rotor-teeter soft-stop linear-spring constant (N-m/rad) [used only for 2 blades and when TeetMod=1]
          0   TeetHSSp    - Rotor-teeter hard-stop linear-spring constant (N-m/rad) [used only for 2 blades and when TeetMod=1]
---------------------- DRIVETRAIN ----------------------------------------------
        100   GBoxEff     - Gearbox efficiency (%) [Envision modification: Alternatively, the name of the file containing gearbox efficiency table]
         97   GBRatio     - Gearbox ratio (-)
8.67637E+08   DTTorSpr    - Drivetrain torsional spring (N-m/rad)
  6.215E+06   DTTorDmp    - Drivetrain torsional damper (N-m/(rad/s))
---------------------- FURLING -------------------------------------------------
False         Furling     - Read in additional model properties for furling turbine (flag) [must currently be FALSE)
"unused"      FurlFile    - Name of file containing furling properties (quoted string) [unused when Furling=False]
---------------------- TOWER ---------------------------------------------------
         20   TwrNodes    - Number of tower nodes used for analysis (-)
"NRELOffshrBsline5MW_Onshore_ElastoDyn_Tower.dat"    TwrFile     - Name of file containing tower properties (quoted string)
---------------------- OUTPUT --------------------------------------------------
True          SumPrint    - Print summary data to "<RootName>.sum" (flag)
          1   OutFile     - Switch to determine where output will be placed: {1: in module output file only; 2: in glue code output file only; 3: both} (currently unused)
True          TabDelim    - Use tab delimiters in text tabular output file? (flag) (currently unused)
"ES10.3E2"    OutFmt      - Format used for text tabular output (except time).  Resulting field should be 10 characters. (quoted string) (currently unused)
          0   TStart      - Time to begin tabular output (s) (currently unused)
          1   DecFact     - Decimation factor for tabular output {1: output every time step} (-) (currently unused)
          0   NTwGages    - Number of tower nodes that have strain gages for output [0 to 9] (-)
         10,         19,         28    TwrGagNd    - List of tower nodes that have strain gages [1 to TwrNodes] (-) [unused if NTwGages=0]
          0   NBlGages    - Number of blade nodes that have strain gages for output [0 to 9] (-)
          0   BldGagNd    - List of blade nodes that have strain gages [1 to BldNodes] (-) [unused if NBlGages=0]
              OutList     - The next line(s) contains a list of output parameters.  See OutListParameters.xlsx for a listing of available output channels, (-)
"BldPitch1"               - Blade 1 pitch angle
"Azimuth"                 - Blade 1 azimuth angle
"RotSpeed"                - Low-speed shaft and high-speed shaft speeds
"GenSpeed"                - Low-speed shaft and high-speed shaft speeds
"TTDspFA"                 - Tower fore-aft and side-to-side displacements and top twist
"TTDspSS"                 - Tower fore-aft and side-to-side displacements and top twist
"TTDspTwst"               - Tower fore-aft and side-to-side displacements and top twist
"RotTorq"                 - Rotor torque and low-speed shaft 0- and 90-bending moments at the main bearing
"LSSGagMya"               - Rotor torque and low-speed shaft 0- and 90-bending moments at the main bearing
"LSSGagMza"               - Rotor torque and low-speed shaft 0- and 90-bending moments at the main bearing
"YawBrFxp"                - Fore-aft shear, side-to-side shear, and vertical forces at the top of the tower (not rotating with nacelle yaw)
"YawBrFyp"                - Fore-aft shear, side-to-side shear, and vertical forces at the top of the tower (not rotating with nacelle yaw)
"YawBrFzp"                - Fore-aft shear, side-to-side shear, and vertical forces at the top of the tower (not rotating with nacelle yaw)
"YawBrMxp"                - Side-to-side bending, fore-aft bending, and yaw moments at the top of the tower (not rotating with nacelle yaw)
"YawBrMyp"                - Side-to-side bending, fore-aft bending, and yaw moments at the top of the tower (not rotating with nacelle yaw)
"YawBrMzp"                - Side-to-side bending, fore-aft bending, and yaw moments at the top of the tower (not rotating with nacelle yaw)
"TwrBsFxt"                - Fore-aft shear, side-to-side shear, and vertical forces at the base of the tower (mudline)
"TwrBsFyt"                - Fore-aft shear, side-to-side shear, and vertical forces at the base of the tower (mudline)
"TwrBsFzt"                - Fore-aft shear, side-to-side shear, and vertical forces at the base of the tower (mudline)
"TwrBsMxt"                - Side-to-side bending, fore-aft bending, and yaw moments at the base of the tower (mudline)
"TwrBsMyt"                - Side-to-side bending, fore-aft bending, and yaw moments at the base of the tower (mudline)
"TwrBsMzt"                - Side-to-side bending, fore-aft bending, and yaw moments at the base of the tower (mudline)
END of input file (the word "END" must appear in the first 3 columns of this last OutList line)
---------------------------------------------------------------------------------------

The problem is that when only the GenDOF is activated the GenIner cannot be set to very small values, and the impact of the GenIner on the rotation behavior of the beam cannot be excluded. Is it possible to let the beam free spinning without the effect of the generator inertia?

Best regards,

[LaurenceWETI]

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,

[jjonkman]

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,

[LaurenceWETI]

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,

[jjonkman]

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,

[LaurenceWETI]

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

[andrew-platt]

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

[LaurenceWETI]

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

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

2. Is there any other method to derive the Eta values?

3. Can you see any mistakes in my calculation?

Best regards,

[LaurenceWETI]

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,

[jjonkman]

I can't spot the error.

[andrew-platt]

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

[LaurenceWETI]

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.

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

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

Best regards

[jjonkman]

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,

[LaurenceWETI]

Dear @jjonkman,

Thank you for your reply.

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

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

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

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

Best regards

[jjonkman]

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,

[LaurenceWETI]

Dear @jjonkman

Thank you for your help

By tracking the values of m%qp%uuu (4:6,:,:) and p%uu0(4:6,:,:), I found that the values of p%uu0 are always the same for blades 1, 2, and 3.

• Should the initial position p%uu0() has the same value for all blades?

• Could that explain, why the transformation is correct for blade 1 and is composed for blades 2 and 3?

Below I run the simulation with only 2 blades.

It is obvious that the transformation for blade 2 uses the same initial position of blade 1. But I don't know how to updatep%uu0()to the value of the second or the third blades.

Best regards,

[andrew-platt]

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