I am trying to make fault affecting blades such as crack, torsion, or bending.
I know for continual deformation I can do it changing input files. Now, I wonder if- and how- I am able to make such changes
during simulation as an abrupt fault; I think it is possible through changing Subroutine Codes provided in FORTRAN, but, actually,
I am not sure, or I think maybe there are some other alternative.

I would be very much obliged if you can help me.
Best Regards

Without changes to the source code, FAST cannot model a sudden blade failure through cracking, delamination, etc., which I presume you’d want to model as a sudden change in blade stiffness. I have not attempted this myself, but one challenge I foresee will be that such blade properties are defined as parameters (constants) at the beginning of the simulation that are not meant to be changed. You can probably get around this by adding a scale factor that changes the stiffness at a given time in the simulation. Perhaps a bigger problem is that the blade mode shapes are derived from the stiffness (among other factors), so, changing the stiffness should also impact the mode shape, which will be more difficult to adjust in the middle of a simulation.

I also want to realize the sudden blade failure, such as changes of the blade stiffness, and the corresponding mode shape of the blade in the middle of the FAST simulation, as you have mentioned before. Now I found the variables in the elastoDyn.f90., i.e. FlpStff,EdgStff, GJStff,(these variables are input by “NRELOffshrBsline5MW_Blade.dat” for linear interpolation, if my understanding is right) The first thing is to change the values of blade properties and blade mode shapes in “NRELOffshrBsline5MW_Blade.dat”. I have four questions.
(1) You mentioned the mode shape is influenced by the stiffness, so how to calculate/determine the values of the mode shape and modify the original values in the “NRELOffshrBsline5MW_Blade.dat” according to the stiffness of blades. I also don’t know what’s the meaning of “BldFl1Sh(2)”, “coeff of x^2”… The original mode shapes are presented as follows.
---------------------- BLADE MODE SHAPES ---------------------------------------
0.0622 BldFl1Sh(2) - Flap mode 1, coeff of x^2
1.7254 BldFl1Sh(3) - , coeff of x^3
-3.2452 BldFl1Sh(4) - , coeff of x^4
4.7131 BldFl1Sh(5) - , coeff of x^5
-2.2555 BldFl1Sh(6) - , coeff of x^6
…
(2) How to consider the stiffness changes appear in the middle of simultion (for example, the blade fails at 10s simulation time)? Can we only modify the source code (elastoDyn.f90)? for instance, add an if-statement to judge the simulation time when stiffness changes appear？Do you have any suggestions?

(3) You mentioned that it is challenging to adjust the mode shape in the middle of simulation because it would be influenced by the stiffness. Likewise, can we also modify the elastoDyn.f90 and add an if-statement to determine the changed mode shape? If it is not possible, could you give me some advice?

(4) If it is difficult to consider the changed mode shape in the middle simulation, I am also wondering if it is possible to divide the one simulation case into two cases, instead of directly consider the stiffness changes in one simulation case? The first case only simulates the normal operation of wind turbine and end up with a fault (such as end up at 10s simulation time). Then the initial boundary condition of the second case is derived from the end of the first case (at 10s simulation time) and begin with the same fault as the first case.

Yes, the mode shapes depend on the distributed mass and stiffness and the appropriate boundary conditions. To derive the blade mode shapes, we normally consider the blade cantilevered to a hub with rigid support structure and derive the mode shapes for a specific rotor speed and blade-pitch angle. BldFl1Sh(2-5) (etc. for each flapwise and edgewise mode) are coefficients of a 6th-order polynomial (without the offset and linear term) that define the mode shape. The NREL preprocessor Modes can calculate the polynomial coefficients directly. The NREL preprocessor BModes can calculate the mode shapes, but you must fit a polynomial to the mode shape to derive the polynomial coefficients (e.g. using the ModeShapePolyFitting.xls spreadsheet). These topics have been discussed many times on this forum.

2/3) Yes, you would have to change the ElastoDyn source code if you want the stiffness and/or mode shapes to change within a given simulation. The input values are used to define variables that are currently stored as time-invariant parameters in the source code. You would have to change these so that the variables are time-dependent parameters. I would start be seeing how the stiffness and mode shapes specified in the ElastoDyn input file are processed by the source code. That said, I’m not really sure why you’d want to to model the actual transient event of a stiffness change. Normally a tool like FAST is used to calculate structural loads that could help you predict the onset of failure (or to design the blade so that it does not fail), but not to directly model an actual blade failure event.

Yes, you could always change the blade stiffness and mode shapes between separate FAST simulations, as long as you don’t want to model the actual transient event associated with the stiffness change.

Thanks for your answers. I am still trying to understand the meaning of the source code in ElastoDyn.f90 and want to make the stiffness or/and mode shapes of the blades changing with time firstly. But now I still confused, I am not sure which part or which subroutine I should modify. Could you give me some suggestions about the specific subroutine and clues? Many thanks to you. Eclosed please find the public subroutines contained in the ElastoDyn.f90.

! … Public Subroutines …

PUBLIC :: ED_Init ! Initialization routine
PUBLIC :: ED_End ! Ending routine (includes clean up)

PUBLIC :: ED_UpdateStates ! Loose coupling routine for solving for constraint states, integrating
! continuous states, and updating discrete states
PUBLIC :: ED_CalcOutput ! Routine for computing outputs

PUBLIC :: ED_CalcConstrStateResidual ! Tight coupling routine for returning the constraint state residual
PUBLIC :: ED_CalcContStateDeriv ! Tight coupling routine for computing derivatives of continuous states
PUBLIC :: ED_UpdateDiscState ! Tight coupling routine for updating discrete states

PUBLIC :: ED_JacobianPInput ! Routine to compute the Jacobians of the output (Y), continuous- (X), discrete-
! (Xd), and constraint-state (Z) equations all with respect to the inputs (u)
PUBLIC :: ED_JacobianPContState ! Routine to compute the Jacobians of the output (Y), continuous- (X), discrete-
! (Xd), and constraint-state (Z) equations all with respect to the continuous
! states (x)
PUBLIC :: ED_JacobianPDiscState ! Routine to compute the Jacobians of the output (Y), continuous- (X), discrete-
! (Xd), and constraint-state (Z) equations all with respect to the discrete
! states (xd)
PUBLIC :: ED_JacobianPConstrState ! Routine to compute the Jacobians of the output (Y), continuous- (X), discrete-
! (Xd), and constraint-state (Z) equations all with respect to the constraint
! states (z)

PUBLIC :: ED_GetOP ! Routine to pack the operating point values (for linearization) into arrays

SUBROUTINE Coeff() in ElastoDyn.f90 is used to set the generalized elastic stiffness matrices (p%KBF for flapwise and p%KBE for edgewise) based on the specified distributed stiffness and mode shapes. You’ll want to change these stiffness parameters (and perhaps mode shapes) to be time dependent, so, you’ll have to change how they are defined/stored (to add a dimension to the array indicating the time variation) and used (to make use of the extra dimension).

With your help, I have added an additional dimension of the variable KBF, KBE and make it changed with time. I am trying to understand the infuence of the changed stiffness and mode shape. According to the ElastoDyn, the KBF and KBE are used to calculate the natural frequency of the blade FreqBE (Inv2Pi * sqrt (stiffness/mass)), and then to calculate the damping CBF. The variable CBF is finally employed in the subroutine FillAugMat to calculate the mass matrix AugMat (if it is right). My question is:
According to the source code, I am wondering whether such a failure (reduced stiffness and corresponding changed mode shape) will only influence the mass matrix and finally influence the motions of the wind turbine. (2) How would the turbine behave if the stiffness is changed? I am not sure how to know my modification of the source code is correct.

Actually, it is still not clear for me to understand the transfer process of the influence of the reduced stiffness, and ask for your help. Thanks for your kind help and timely responses all the time!

KBF and CBF are the generalized elastic stiffness and structural damping matrices, respectively. They do not influence the mass matrix, but rather they influence the forces on the right-hand side of generalized equations of motion, m*a=F (AugMat contains the generalized mass matrix and force vector concatenated into a single matrix, [m F]).

If you are changing KBF (which then changes CBF) during the course of a simulation, I would expect the natural frequencies and deflections for a given load to change. You could do a sanity check if this is working by:
(1) disabling all features of FAST except ElastoDyn,
(2) disable all DOFs of ElastoDyn except the blades,
(3) set the initial rotor speed to zero,
(4) set a nonzero initial blade displacement (nonzero OoPDefl), and
(5) simulate this free-decay motion.

I would expect a noticeable change in oscillation frequency and motion amplitude of the blade deflection outputs after changing KBF (likewise for KBE).

Many thanks to you for your kind help all the time. In addition to the blade failure, I also want to make the mooring line failure. I am considering two outmost failures now. One is a failure in the top segment (fairlead), which is close to the floating platform (hence the mooring line will fall away and the tension will not pull anymore the floater if it happens), another one is a failure in the bottom segment of the line which is close to the anchor (the anchor will move and there is a larger part of the line resting on the seabed).

Similarly, I hope to revise the code of the MoorDyn to make the fault. In my idea, (1) I thinkg reducing the tension of the top segment of the line into a small value or zero, is enough to represent the first failure. But I am not sure. (2) As for the second failure, is it right if I directly change the position of the bottom node of the mooring line in the MoorDyn code?

Simulating mooring line failures would be an interesting capability. I’ve thought of the idea itself, but never how to go about realizing it in MoorDyn. Here are a few initial thoughts in case they are of help.

Most of MoorDyn’s information about the mooring system is stored in “m” (the MiscVar structure), so you would probably want to add variables within it to specify the details of the failure, unless you chose to hardcode in the specifics of each failure.

I can think of two general approaches for implementing the failure:
(1) Modifying the data in m once the failure occurs to represent the disconnection without changing MoorDyn’s calculation process, or
(2) Modifying the calculations within MD_CalcContStateDeriv (specifically within DoLineRHS) to get the same effect as a failed line.

Your ideas are in the second category, which is I think a good place to start. For the fairlead failure, I agree that zeroing the tension seems reasonable. You could simply set FairFtot to zero in DoLineRHS for the failed line. For the anchor failure, you could change the anchor node position (Line%r(J,0) in DoLineRHS), but abrupt position changes might cause some problematic damping forces. An option that I think would be just as effective and potentially simpler would be to set the length of the line segment that attaches to the anchor to be very large (i.e. set Line%l(1) in DoLinesRHS to a large number).

Hope that helps, and I look forward to hearing how it goes!

Thanks for your kind suggestions, which really help me a lot. I am sorry for not response to you so quickly, because I am working on the problem of the mooring line failure these days. Following your advice, I set the m%LineList(1)%FairConnect)%Ftot to be zero and set the m%LineList(1)%l(1) to be a large value (e.g. 3000m) in the subroutine of MD_CalcContStateDeriv in front of the subroutine DoLineRHS (I guess modifying the code in the DoLineRHS is also fine, but I am not sure if it could consider the failure of one/two mooring line(s)…), which is the simplest method you mentioned. It seems it workes. Attached please find some preliminary results of the mooring line failure that I want to share with you. The failure time is 100s. Many thanks to you.

I am trying to add some band-limited white noises into the output signals (e.g. tower top acceleration y) to represent the measurement noise of tower top acceleration y sensor. The magnitude of the noise is determined according to the following literature (Table 1) which has been enclosed: Odgaard, P., Johnson, K., Wind turbine fault detection and fault tolerant control - a second challenge.
According to the paper, the noise power of the band-limited white noise for the tower top acceleration y is 5*10^-4. It means that the magnitude of the noise is around ~0.1- 0.15m/s^ 2 if the sampling time is 0.1 according to the band-limited white noise block. But I found the tower top acceleration y derived from FAST simulation is quite small, only around 0.02-0.04m/s^ 2 (attached). it would be impossible to detect the signal of tower top acceleration y if I add such a noise. Since I don’ t know the real values of the measurement noise of the sensor, I am not sure how to consider the problem and ask for your help. Can you help me? Thank you.

In terms of configuration of ElastoDyn module, I only activate the First tower bending mode DOF for simplification as follows. The second tower bending mode DOF is disabled in order to develop another simple linear model for offshore wind turbine.
True FlapDOF1 - First flapwise blade mode DOF (flag) !
True FlapDOF2 - Second flapwise blade mode DOF (flag) !
True 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) !
True TwFADOF1 - First fore-aft tower bending-mode DOF (flag)
False TwFADOF2 - Second fore-aft tower bending-mode DOF (flag) !
True TwSSDOF1 - First side-to-side tower bending-mode DOF (flag)
False TwSSDOF2 - Second side-to-side tower bending-mode DOF (flag) !
True PtfmSgDOF - Platform horizontal surge translation DOF (flag)
True PtfmSwDOF - Platform horizontal sway translation DOF (flag)
True PtfmHvDOF - Platform vertical heave translation DOF (flag)
True PtfmRDOF - Platform roll tilt rotation DOF (flag)
True PtfmPDOF - Platform pitch tilt rotation DOF (flag)
True PtfmYDOF - Platform yaw rotation DOF (flag) 10.1109_ACC.2013.6580525.pdf (370 KB)

I can’t really comment on the magnitude of expected measurement noise.

Without knowing more about your simulation setup, it is hard to know why your tower-top acceleration in the side-to-side (y-direction) is so low. What wind turbine are you modeling? Are you consider normal operation under turbulence? Are you getting the side-to-side deflections and velocities that you’d expect?

Thank you for your quick response. I use a DTU 10MW reference wind turbine and the floating foundation is SWE-TripleSpar Floating
Platform (attached). It is in normal operating condition with constant wind speed (12m/s) without turbulence, irregular wave based on JONSWAP spectrum with height 1.26m and period 10s in the streamwise direction. The natural frequency of the whole turbine is:
Mode Natural frequency [Hz]
1st Tower fore-aft mode 0.25
1st Tower side-side mode 0.25
1st collective flap mode 0.63

I’m not too familiar with the DTU 10-MW model, but the steady-state tower-top deflections look reasonable to me. It looks a bit odd to me that the tower side-to-side doesn’t start deflecting until about 50 s into the simulation. Do you know what in your simulation set-up is causing this?

Regardless, I would not expect large tower side-to-side accelerations without wind turbulence; you should see much higher accelerations when turbulence is included.

It’s very kind of you. Thank you for helping check the results. In my case, the OWT starts to generate wind power only after 50s since the rotor speed and generator speed have a set-up process and increase from 0 at the beginning (0s) and reach rated value at around 50s. I think this is the reason that the tower side-to-side don’t deflect until ~50s. Anyway, good to know that such a steady tower top acceleration seems reasonable and I will try to include some turbulence if large tower acceleration is expected. Thanks again.

Thank you for your help last time and I also hope to realize the same mooring line failure (fairlead and anchor failures) in the MoorDyn C, namely the C lib. that I downloaded on your website (matt-hall.ca/moordyn.html). It also would be connect to the MATLAB in my case. The fairlead failure is zeroing the tension at the fairlead while the anchor failure is realized by a large length of the line segment attached to the seabed, as you suggested. But it seems that the source code in MoorDyn C is a much different with the MoorDyn in FASY v8 that I asked before. I have no idea which source code (Connection.cpp, Line.cpp, NoorDyn.cpp), which variable and which sub-function I can revise. Can you give me some clues? Many thanks to you and look forward to hearing from you.

I’m intending to do similar thing what Yichau.Li already tried to do. I’m trying to change the blade mass properties dynamically within a given simulation. I know such changes is not possible without changing the source code. But I’m just wondering which source code must be changed, to consider such changes, is that the ElstoDyn.f90 or the BeamDyn.f90. Because both of this modules are dealing with blade input file to calculate blade motion and reaction loads. I started to read the BeamDyn.f90 source code to understand where and how the mass matrix are calculated, but till now I did not make any progress. Do you have any advice how this should be implemented

While both ElastoDyn and BeamDyn can be used in FAST / OpenFAST to model blade structural dynamics, the choice of which module is used is selected in the primary input file (CompElast = 1 is for ElastoDyn; CompElast = 2 is for BeamDyn).

I provided some guidance above about how to go about implementing this feature in ElastoDyn in the forum posts in this topic above.

I tried to modify the Blade element mass during the simulation by adding a time condition in the SOUBROUTINE Coeff() as below:

IF (t > 5) THEN ! Add mass to blade Element at simulation time 5 [s]
p%BElmntMass(J,K) = 1.5*p%MassB(K,J)*p%DRNodes(J) ! Mass of blade element J Multiples by 1.5
ELSE
! Calculate the mass of the current element
p%BElmntMass(J,K) = p%MassB(K,J)*p%DRNodes(J) ! Mass of blade element J
ENDIF
Therefore I added a time variable (t) to the SUBROUTINE Coeff() and to all Subroutines, where this Subroutine is called.
But the compiling of OpenFAST after this modification didn’t succeed. I think the problem is that by adding the time variable (t) to the Subroutines, which are only one time called at the beginning of the simulation.
I’m trying now to add a separate Subroutine, where I can add some mass to blade element, but I think I’ll have same problem as above, because I need this parameters to be displayed in the SUBROUTINE SetOtherParameters ()
Do you any advice how this should be implemented? Any help will be really appreciated!