Implementation of external forces in BeamDyn

[LaurenceWETI]

Hello everyone

@jjonkman @bjonkman

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

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

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

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

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

I would appreciate any help.

Best regards

Laurence

[jjonkman]

Hi @LaurenceWETI,

Overall the approach sounds OK. Just a few comments:

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

Best regards,

[LaurenceWETI]

Hello @jjonkman,

Thank you for the quick reply.

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

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

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

Best regards Laurence

[andrew-platt]

Hi @jjonkman and @LaurenceWETI,

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

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

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

Regards, Andy

[jjonkman]

Hi @andrew-platt and @LaurenceWETI,

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

Best regards,

[LaurenceWETI]

Hello @andrew-platt and @jjonkman,

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

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

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

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

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

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

Best regards, Laurence

[andrew-platt]

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

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

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

Regards, Andy

[jjonkman]

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

Best regards,

[LaurenceWETI]

Hello @andrew-platt and @jjonkman,

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

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

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

Best regards,

[andrew-platt]

A few thoughts on this:

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

Regards, Andy

[LaurenceWETI]

Hello @andrew-platt,

Thank you for your support.

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

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

Best regards, Laurence

[LaurenceWETI]

Hello everyone,

@jjonkman @andrew-platt @bjonkman

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

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

This equation is implanted in BD source code to compute the x, y and z components of the distributed Coriolis force p%DistrCoriolisForce at node number 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,

[jjonkman]

Dear @LaurenceWETI,

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

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

Best regards,