Hi,
I am attempting to use MLIfe with a user defined, 4 variable probability distribution and I have ran into some questions/problems.
When selecting to use a User Defined Distribution, why is it necessary to still specify the Wind Speed Channel for fatigue calculations? Is there any way around this? Why is this necessary if using a User Defined Distribution? I have ran a large number of simulations without Wind Speed as an output and would ideally avoid re-running them!
Thanks,
Luke Cunningham
Dear Luke,
The mean wind speed of each simulation is needed in MLife’s fatigue calculation to determine which wind-speed bin the simulation data falls in to assign the proper probability of that simulation for its contribution to the overall lifetime.
If your simulation data does not have a wind speed output, presumably you know in some other way which bin each simulation is associated with. But you’d have to modify the MLife source code to define this without the wind speed output.
Best regards,
Hi Jason,
Thank you for your reply. I now have Mlife running with the User Defined Distrubtion.
Perhaps you could offer further insight into some questions I have regarding using Mlife:
When defining th LUlt for parameters, is it assumed that the LUlt will be in the same units as the paramter is defined? For example, FAST outputs 'FAIRTEN’in N, while a parameter such as ‘TwrBsFxt’ is listed in KN. I assume the LUlt for FAIRTEN should be defined in N and that of TwrBsFxt in KN?
What would be common reason for the ‘Lifetime Damage’ and 'Time until Failure 'to be calculated as negative numbers or ‘NaN’? The time series I have develoepd for the outputs are all fully intact (i.e. no NaN values) and look to me to be reasonable. I would have thought that even if parameters such as ‘LUlt’ or ‘m’ are poorly chosen that this would result in very small or very large values of ‘Lifetime Damage’ and ‘Time until Failure’ and, as such, would assume there is something wrong with my .mlif file inputs. I have attached it below as reference. TCDCoupledExample.txt (109 KB)
Thank you for your time and assistance, as always,
Kind regards,
Luke Cunningham
Dear Luke,
Regarding (1), yes.
Regarding (2), I’m not sure. I don’t recall running into that problem myself (granted, I don’t run MLIfe very often). I would probably start by reviewing the rainflow cycle counts for correctness.
Best regards,
Hi Jason,
Okay, thanks. I will have a look at that.
Thank you,
Luke Cunningham
Hi,
I am conducting further work based on the research Luke carried out. I am having a similar issue to what he describes in the second question in his second post:
Looking at the short term fatigue statistics I have discovered that several of the 404 different cases that I ran feature the value ‘inf’ for the tower base moments in the x and y directions, so I assume this leads to the ‘NaN’ value when extrapolating for the lifetime fatigue statistics. I have looked extensively through the MLife code which computes the short term fatigue statistics and cannot find the reason behind this ‘Inf’ value. I think the ‘inf’ value may be arising because somewhere in the code a number is being divided by 0. However I cannot find the source of it and I am confused as several of my test cases compute the fatigue statistics correctly while other cases do not (each case was simply an OpenFAST simulation with varying wind and wave conditions).
Like Luke I am using the same 4 variable probability distribution however I don’t think this would have an effect on the short term fatigue statistics if it were to be incorrect. I am starting to lean towards the fact that this code be an issue with the MLife code?
I am wondering if anyone has a solution to this problem as I have run out of ideas. I have attached the short term fatigue statistics which show how some cases are not being computed correctly while others are. I would be vary grateful for any help on this issue.
Kind regards,
Conor White.
TMD_Short-term_Damage_Rate.txt (196 KB)
Dear Conor,
What does the time series look like for a case that is resulting in “Inf” in the short-term damage rate? Do you have the same problem whether you enable or disable the Goodman correction?
Best regards,
Hi Jason,
Thank you for the fast response. Only one of the 404 output file produces ‘inf’ for the short term fatigue for the tower base moments in the x direction TwrBsMxt. I have attached this file called ‘OFWT_11.00_114.28_3.95_12.73.out’. Interestingly this is the only output file of the 404 which produces the ‘inf’ value for the TwrBsMxt (each of the 404 output files are just varying metocean conditions).
I have also attach the file “OFWT_20.10_-24.69_3.95_12.73.out” which is an arbitrarily chosen example of one of the cases in which the ‘inf’ value occurs for the tower base moment in the y direction. Of my 404 cases about 100 produce this ‘inf’ value for the TwrBsMyt. All of these time-series appear to be correct. Please note that I have only attached a small portion of the two output files as the original hour long time-series was too large to attach.
The problem doesn’t occur when I disable the Goodman correction so obviously the issue must be arising from the part of the code that computes this correction. I have examined this in the source code however the solution isn’t obvious to me as none of the variables in the ‘cyclesToFailure’ equation appear to produce values that would result in ‘inf’ for the short-term fatigue statistics. I have included the code extract which I am referring to from the ‘compute_fatigue_per_DEL_type.m’ function:
switch DEL_Type
case 1 % Goodman correction with a specified fixed mean
cycleRanges = cycles(:,3);
case 2 % Goodman correction with a fixed zero mean
cycleRanges = cycles(:,3); % was 5
case 3 % no Goodman correction
cycleRanges = cycles(:,1);
otherwise
% error
end
And later in the code:
% Compute cycles to failure per Equation 6 of the Theory manual.
cyclesToFailure = ( ( Fatigue.ChanInfo(iCh).LUlt - abs( double(lmf) ) )./( 0.5*cycleRanges ) ).^SNslope;
UPDATE:
I just realised could this issue be due to choosing an incorrect value for LUlt, given that it is used in the cyclesToFailure equation above? Given I am computing the tower base moments and shear force I simply used the value for LUlt in which Luke (the previous poster - I am conducting similar work to him) used which was 31,875 kN. I haven’t yet verified where he got this value from so it may be incorrect.
Thanks for the help,
Conor.
OFWT_20.10_-24.69_3.95_12.73.txt (60.9 KB)
OFWT_11.00_114.28_3.95_12.73.txt (123 KB)
Dear Conor,
Yes, my guess is the problem is related to the setting of LUlt. It looks like many of the loads are higher than 31875, which would indicate immediate failure. I would recommend setting LUlt above the highest absolute maximum load across all time series (and distinct for shear forces, axial forces, bending moments, and torsion). See our related FAQ for more information: nwtc.nrel.gov/faq#TypeLMF.
Best regards,
Hi Jason,
Thank you for your reply. setting LUlt above the highest absolute maximum load across all time series solved the issue, however I have another question;
Given that I am looking at the tower base shear forces in the x,y and z directions and the tower base moments in the x,y and z directions I used MExtremes to find the absolute maximum values for these 6 channels across my 404 .out files. I applied the methodology of scaling these maxima by a factor of 20 to obtain the values for LUlt as suggested in the FAQ. My question is should I take the maximum value of the three base moment channels x,y and z for example, and use this number multiplied by a factor of 20 for the LUlt for all of the base moment components (i.e TwrBsMxt, TwrBsMyt and TwrBsMzt) or should I take the maximum of each individual base moment channel scaled by a factor of 20 and use that as the LUlt value? e.g using three different LUlt values for the three base moment channels as opposed to using the overall maximum LUlt value for the three base moment channels.
Thanks for the help,
Conor.
Dear Conor,
I’m glad increasing LUlt solved your problem.
Assuming your tower is axisymmetric, I would set 4 unique values of LUlt for each cross section, i.e.:
Tower-Base Shear Forces (x or y): FactorMAX( SQRT( TwrBsFxt^2 + TwrBsFyt^2) )
Tower-Base Axial Force (z): FactorMAX( ABS(TwrBsFzt) )
Tower-Base Bending Moments (x or y): FactorMAX( SQRT( TwrBsMxt^2 + TwrBsMyt^2) )
Tower-Base Torsion Moment (z): FactorMAX( ABS(TwrBsMzt) )
Note that MAX( SQRT( x^2 + y^2) ) is not the same as SQRT( MAX(x^2) + MAX(y^2) ).
Note also that you don’t need to use Factor = 20, but presumably Factor > 1 to avoid immediate failure. If you know the ultimate strength of the cross section, you could derive LUlt from that.
Best regards,
Hi Jason,
I now understand what you mean about setting 4 unique values for LUlt, thank you for the advice.
My research is concerned with the 5MW baseline OC3-Hwyind tower, to my knowledge there is no specified value for the ultimate strength of the tower. However, from other posts I have read that in order to determine the ultimate strength for the bending moment the felxure formula can be used (LUlt = sigma_yield*Ixx/y_max), where y_max was taken to be diameter/2. Given that the diameter of the Hwyind tower is 6.5m, its thickness is 0.027 m and E = 200 GPA (I found the values here nrel.gov/docs/fy10osti/47535.pdf) , I calculated LUlt to be 372,184.62 [kN m].
However when I use this value for LUlt within MLife (for TwrBsMxt and TwrBsMyt) I get a result fatigue life of approximately 9.85e+04 seconds which clearly seems incorrect given that the tower has a 20 year design life. I am quite confused with this result given that my simulated data (within TurbSim and OpenFAST) is correct to my knowledge.
I was wondering does my calculation of LUlt look correct to you and could you suggest any reason for the issues with the lifetime fatigue statistics?
Thanks for the help,
Conor White.
Dear Conor,
I agree with the equation LUlt = sigma_yield*Ixx/y_max for the bending moment. But using your values, I obtain Ixx = 2.876 m^4 and LUlt = 1.77E8 kNm.
Best regards,
Hi Jason,
Apologies I made an error in my calculations and I agree with the value you obtained for LUlt. I noticed that you have used 200 Mpa or 2E+8 N/m^2 as the yield strength for the material - is this specified for the OC3-Hywind tower? I ask because I was under the impression that the yield strength of steel can vary widely (200 - 2100 MPa).
Using this value for LUlt produces more reasonable lifetime statistics for the TwrBsMxt and TwrBsMyt - am correct in assuming this value of LUlt can be used for both these bending moments? The TwrBsMzt on the other hand, I am unsure whether this LUlt value can used?
As for LUlt for the shear force, I have read else where on the forum that simply, sigma = F/A is used. Am I correct in saying that in this case LUlt for the TwrBsFxt, TwrBsFyt an TwrBsFzt can be computed by rearranging this equation to give LUlt = sigma_yield*A ? If this is correct am I right in saying that the area, A would be the area of the cross section of the tower base?
I greatly appreciate your help, it has cleared up a lot for me.
Regards,
Conor.
Dear Conor,
Just to be clear, the OC3-Hywind specifications report (nrel.gov/docs/fy10osti/47535.pdf) gives the Young’s modulus (E) of the tower as 210 GPa; however the yield strength (sigma_yield = E*epsilon_yield, where epsilon_yield is the axial strain at the yield strength) is not given in that report. I misread your post and was simply using your value. You’ll have to make your own assumptions on sigma_yield for your own purposes.
Here are the basic equations for LUlt for a thin-walled circular cross section (from basic stress equations and tau_yield = G*gamma_yield, where G is the shear modulus and gamma_yield is the shear strain at the yield strength; you’ll have to assume your own value of tau_yield):
Transverse shear: Lult = tau_yield * A/2
Axial loading: LUlt = sigma_yield * A
Pure bending: LUlt = sigma_yield * I/y
Torsion: LUlt = tau_yield * J/r
where:
Outer diameter: D_o
Inner diameter: D_i
Maximum radius: y = r = D_o/2
Area: A = pi/4 * ( D_o^2 - D_i^2 )
Area transverse inertia: I = pi/64 * ( D_o^4 - D_i^4 )
Area polar inertia: J = 2 * I = pi/32 * ( D_o^4 - D_i^4 )
I hope that helps.
Best regards,
Dear Jason,
I encountered the same problem as depicted by Conor, and the calculated failure time of the tower base is limited to 3.78E+05s, which is much lower than the specified service time, 20 years (6.3072E+08s). The detailed file is uploaded to my Google cloud drive, considering the large input file with a simulation time of one hour. drive.google.com/drive/folders/ … sp=sharing
I listed the relevant main setting parameters in the “.mlife” file to clarify the problem.
WeibullShape: 1.93
WeibullScale: 5.71
FiltRatio: 0.0
UCMult: 0.5
BinCycles: TRUE with BinFlag=BN and BinWidth/Number=20
SNslopeLst=10 for the blade and SNslopeLst=3 for the tower
TypeLMF=AM
LUlt is the key parameter for the real fatigue lifetime prediction according to my experience, and I adopted the following information to consider the relevant values given that all the component is under the pure bending states even though this assumption is not realistic but applicable.
For the blade,
Because detailed layup information is unavailable, I referred to the 5MW/61.5 blade designed by Sandia National Laboratory (osti.gov/biblio/1095962). And the detailed information and corresponding results are shown below:
Outer diameter 3.542 m
thickness 0.05457 m
Inner diameter 3.43286 m
Moment of Inertia 0.909152339 m^4
Tensile strength 700000 kPa
Ultimate bearing capacity 359348.7504 kN-m
For the tower, I adopted the parameters in the reporter entitled “Definition of a 5-MW reference wind turbine for offshore system development”
Outer diameter 6 m
thickness 0.027 m
Inner diameter 5.946 m
Moment of Inertia 2.259488151 m^4
Tensile strength 470000 kPa
Ultimate bearing capacity 353986.477 kN-m
As I did not find any mistake, and pls help me check the file if possible. Many thanks.
Best regards,
Jian
Dear Jian,
Overall your process looks correct to me. One detail you missed was the tower thickness scaling factor of 1.3, so the tower-base thickness should actually be 0.027*1.3 = 0.0351 m, which will increase LUlt of the tower-base to 458319 kNm. What does that do to the tower design life?
Best regards,
Dear Jason,
Thank you for your reply. I changed the ultimate loading bearing capacity to 458319 kN-m. The failure time increased from 3.78E+05s to 8.84E+05s (around 0.03 years) but did not satisfy the design demand. I first thought it might be caused by the parameter setting in Turbsim, which needs to be improved but is not finished regarding decouple hub height from the grid. (github.com/OpenFAST/openfast/issues/199)
The following part demonstrates the detailed setting in TurbSim (6m/s), which is employed to generate a 360m-height grid employed in FAST.Farm
--------Turbine/Model Specifications-----------------------
37 NumGrid_Z - Vertical grid-point matrix dimension
201 NumGrid_Y - Horizontal grid-point matrix dimension
0.05 TimeStep - Time step [seconds]
4000 AnalysisTime - Length of analysis time series [seconds]
“ALL” UsableTime - Usable length of output time series [seconds] (program will add GridWidth/MeanHHWS seconds) [bjj: was 630]
185 HubHt - Hub height [m] (should be > 0.5GridHeight)
360 GridHeight - Grid height [m]
2000 GridWidth - Grid width [m] (should be >= 2(RotorRadius+ShaftLength))
0 VFlowAng - Vertical mean flow (uptilt) angle [degrees]
0 HFlowAng - Horizontal mean flow (skew) angle [degrees]
--------Meteorological Boundary Conditions-------------------
IECKAI TurbModel - Turbulence model (“IECKAI”=Kaimal, “IECVKM”=von Karman, “GP_LLJ”, “NWTCUP”, “SMOOTH”, “WF_UPW”, “WF_07D”, “WF_14D”, or “NONE”)
“unused” UserFile - Name of the file that contains inputs for user-defined spectra or time series inputs (used only for “USRINP” and “TIMESR” models)
“1-Ed3” IECstandard - Number of IEC 61400-x standard (x=1,2, or 3 with optional 61400-1 edition number (i.e. “1-Ed2”) )
14 IECturbc - IEC turbulence characteristic (“A”, “B”, “C” or the turbulence intensity in percent) (“KHTEST” option with NWTCUP, not used for other models)
NTM IEC_WindType - IEC turbulence type (“NTM”=normal, “xETM”=extreme turbulence, “xEWM1”=extreme 1-year wind, “xEWM50”=extreme 50-year wind, where x=wind turbine class 1, 2, or 3)
default ETMc - IEC Extreme turbulence model “c” parameter [m/s]
“PL” WindProfileType - Wind profile type (“JET”=Low-level jet,“LOG”=Logarithmic,“PL”=Power law, or “default”, or “USR”=User-defined)
“unused” ProfileFile - Name of the file that contains input profiles for WindProfileType=“USR” and/or TurbModel=“USRVKM” [-]
90.0 RefHt - Height of the reference wind speed [m]
6.0 URef - Mean (total) wind speed at the reference height [m/s]
default ZJetMax - Jet height [m] (used only for JET wind profile, valid 70-490 m)
default PLExp - Power law exponent [-] (or “default”)
default Z0 - Surface roughness length [m] (or “default”)
Even though the abovementioned setting is not precise, I checked the maximum bending moment at the tower base between cut-in wind speed and cut-out wind speed, around 100.4 MN-m, which is a rational value to my knowledge. Then I have no idea about the reason causing such a short failure time 
. Do you have any other suggestions? Many thanks.
Best regards,
Jian
Dear Jian,
Your TurbSim inputs look reasonable to me. Presumably you’ve also set the reasonable spatial coherence parameters, as discussed in section 5.2.3 of the FAST.Farm User’s Guide and Theory Manual: nrel.gov/docs/fy21osti/78485.pdf.
I’m not really sure why the tower fatigue lifetime you are predicting is so low. The fatigue lifetime of tower of the NREL 5-MW baseline turbine was never checked when the baseline model was established (the tower was more prescribed based on industry trends rather than designed), but I don’t think the tower details are way off the mark.
Best regards,
Dear Jason,
Many thanks. I consider the coherence parameters (specified as “default”) might be right when not considering wake meandering, as shown below.
--------Spatial Coherence Parameters----------------------------
default SCMod1 - u-component coherence model (“GENERAL”, “IEC”, “API”, “NONE”, or “default”)
default SCMod2 - v-component coherence model (“GENERAL”, “IEC”, “NONE”, or “default”)
default SCMod3 - w-component coherence model (“GENERAL”, “IEC”, “NONE”, or “default”)
default InCDec1 - u-component coherence parameters for general or IEC models [-, m^-1] (e.g. “10.0 0.3e-3” in quotes) (or “default”)
default InCDec2 - v-component coherence parameters for general or IEC models [-, m^-1] (e.g. “10.0 0.3e-3” in quotes) (or “default”)
default InCDec3 - w-component coherence parameters for general or IEC models [-, m^-1] (e.g. “10.0 0.3e-3” in quotes) (or “default”)
default CohExp - Coherence exponent for general model [-] (or “default”)
I am also wondering why the failure time is so low. I tried another method adopting a factor of 43, which corresponds to the ultimate bearing capacity of 5,306,200 kN-m (even though the designed bearing capacity is 458,319kNm based on the material and section properties) and lifetime is 52.6 years. It sounds rational. I guess whether some additional stiffener to the tower base could cause larger bearing capacity in the real design process—just my guess. I will check other parameters to find the problem.
Thank you very much.
Best regards,
Jian