Turbine-soil interaction. Influence on the mode shapes.

Dear Jason,
Thanks for your explanation. it’s my mistake to mix BModes and Modes… :blush:

Dear Sail,

FYI: We now recommend that BModes be used in place of Modes. BModes employs a higher fidelity model than Modes and also includes the option for tower-top inertia (in addition to point mass) and tower-base DOFs for foundation flexibility or floating platforms. The biggest issue with BModes is it doesn’t directly output the polynomial coefficients needed by FAST. We’ve supplied the “ModeShapePolyFitting.xls” spreadsheet in the FAST archive to aid in this effort of deriving the polynomial coefficients.

Best regards,

Dear Jason,
I am modelling the 1.5-MW 3-bladed WT (described in the WindPACT Turbine Rotor Design Study) with my finite element-boundary element code.
I need to know the dimensions of the tapered tower.
Reading the report of the “WindPACT Turbine Rotor Design Study” I found the following geometry from the “Table 6.2-Summary of Properties of the Baseline and Task #5 Final Configurations” :

Thick_top= 0.0087 m
Thick_base= 0.0174 m
D_base= 5.663 m
D_top= 2.823 m
while from the Test13.fst file I found out that:
Flexible Heigth of the tower= H= 82.39 m
I assumed the following steel properties:
Ex= 2.00E+11 N/m^2
Ey= 2.00E+11 N/m^2
G= 7.70E+10 N/m^2
rho= 7.85E+03 kg/m^3
Poisson´s ratio= 3.00E-01
I assumed at first that for top and base diameters the author meant the internal diameters of the top and base annulus and I calculated the distributed properties of the tower (included in the attachment).
The volume of the tower is H*(A_top+A_base)/2=1.60E+01 m^3 while the tower mass is Volume*rho=1.25E+05 kg
My results are in disagreement with the ones recorder in Baseline_Tower.dat.
Even if I assume that the top and base diameters are meant as external diameters, the results do not match the example in The CertTest folder.
However I decide to substitute the values in the file Baseline_Tower.dat with my results and when I run the program I obtain a funky value of the tower mass:
(From the fsm file)
Tower Mass (kg) 116393.037
which is different than 1.25E+05 kg.
Am I doing something wrong? Maybe I confound myself but I do not see the mistake in my calculations.
Are my assumptions for the material properties right? Maybe the density of the steel is assumed different than rho=7.85E+03 kg/m^3.

Thank you for any help.
Best Regards,

Dipl. Ing. Francesca Taddei

RWTH-Aachen University
Lehrstuhl für Baustatik und Baudynamik
Tel: +49 (0)241 / 80 25089

Tower properties.txt (2.49 KB)

Dear Francesca,

I wasn’t the one who originally derived the distributed tower properties of the WindPACT 1.5-MW baseline wind turbine supplied with the FAST archive. Instead, these came from data supplied by GEC (one of the leads on the WindPACT studies). I quickly checked the source of the distributed tower properties and see that they were derived from slightly different values of the diameter/thickness than you report. Here are the values that were used:

base diam = 5.61 m
top diam = 2.82 m
base thickness = 0.01756 m
top thickness = 0.00882 m

The other values you report, including the flexible tower height and steel properties are the correct, except that the density was increased by 5% to account for paint, bolts, flanges, etc. that are not accounted for in the thickness data.

Table 6-2 in the “WindPACT Turbine Rotor Design Study” refers to configuration 1.5A08C01V03cAdm. The version used to derive the WindPACT 1.5-MW model supplied with the FAST archive was configuration 1.5A08V07adm, which is a newer version according to the configuration ID described in Table 2-4 of the same report. I’m not sure why the report lists properties of an older version.

I haven’t checked all of the equations in your attachment, but hopefully the changes in tower values described above will enable you to reproduce the distributed tower properties supplied with the FAST archive.

Best regards,

Hello Jason,

I am looking at putting monopiles or other fixed foundation under 13.2MW turbine in GoM shallow waters. You said FAST can do “coupled springs” or “distributed springs” representations of the foundation. I am not sure I got these functions of FAST so far. And how to use that?

Another question is, up to now HydroDyn can use Linear Wave Theory to model waves. For shallow waters, is there any mathematical model existing in HydroDyn such as Stream Function Theory for calculation.

Thanks so much.

Regards, Lingling

Dear Lingling,

It is not possible in the current version of FAST to directly enable a coupled springs or distributed springs model solely through settings in the input file. Instead, a slight customization (requiring recompile) of FAST is currently required. The coupled springs model can be implemented in FAST through the UserPtfmLd routine, which allows one to apply a discrete load at the platform reference point. The distributed springs model can be implemented in FAST through the UserTwrLd routine, which allows one to apply distributed loads (loads per unit length) along the flexible tower. Please find the UserPtfmLd and UserTwrLd routines that I created for modeling the coupled springs and distributed springs models in the IEA Wind Task 23 Subtask 2 Offshore Code Comparison Collaboration (OC3) project attached. You should be able to use/adapt these attachments to introduce soil flexibility into your FAST model.

The current version of HydroDyn for monopiles accounts for regular or irregular linear waves (with or without stretching) and sea currents and uses the relative form of Morison’s equation for the load calculation. Because the module does not currently have a higher-order wave kinematics model built into it (e.g., Stream Function theory), if you want to model severe regular waves, you would have to use GH Bladed or some other wave kinematics code to generate the higher-order wave kinematics data before running a simulation with HydroDyn. The WaveMod = 4 feature was added so that HydroDyn could read data generated externally by GH Bladed or an equivalent code.

Best regards,
UserTwrLd_DS.f90.txt (10.6 KB)
UserPtfmLd_CS.f90.txt (6.35 KB)

Got you! Thanks!

Dear Jason,

I’m currently working with UserPtfmLd in UserSubs to implement mudline flexibility and damping in FAST for the offshore wind turbine (supported by a monopile). Can you explain what role the platform added mass (PtfmAM) plays in the mudline calculations, and what it defines in a fixed-bottom context? I found that without specifying a rather large platform added mass, my simulations won’t run - and depending on the mass, it adds a “second” first natural frequency to my structure (e.g. a second peak in a Fourier Transform of free vibration aside from the first natural frequency of ~0.3 Hz).

Thank you,

– Wystan

Dear Wystan,

The platform added mass (PtfmAM) within UserPtfmLd() could be used to represent acceleration-dependent loads transfered to the pile from the soil. However, I’ve seen soil models that ignore these terms; so, you can probably set PtfmAM to zero.

Enabling the platform degrees-of-freedom (DOFs) will introduce additional natural modes in the FAST model. My guess is the model is not running because these new modes are at a high enough frequency that that the structural time step must be reduced to resolve them without numerical instability. (And introducing PtfmAM lowers the natural frequency of these modes, making the model numerical stable.) See my post dated Nov 4, 2010 in the forum topic found here for rules of thumb on how to choose proper time steps: http://forums.nrel.gov/t/error-in-fast-working-in-adams/298/1.

Best regards,

Dear all,

a short question: If I add soil matrices (stiffness and damping) in the user-subroutine, in the case of soft soil, FAST delivers very good and reasonable results. If I implement very stiff soil, which means very high stiffness coefficients, the program aborts saying…
I though about a convergence problem and I reduced the time step, succeeding with my stiff soil without problem.
Therefore my conclusion is that, if I add a very stiff soil and with it a high stiffness, the algorithm starts to have problem with convergence as the system is characterized somehow by high frequencies (and low periods) and the time step need to be reduced.
Does it make sense? Is that the physical reason?


Dear Francesca,

Your understanding is correct.

See my post dated Nov 4, 2010 in the forum topic found here for rules of thumb on how to choose proper time steps: http://forums.nrel.gov/t/error-in-fast-working-in-adams/298/1.

Best regards,

Thank you Jason,
now it is clear.

Dear all,

If I want to incorporate p-y curve in the UserTwrLd subroutine, do I have to specify the resistance (p) equation related to displacement (y) just right at the TwrNode position where the displacement vector is calculated?

Did somebody have such examples to share with me? I don’t know whether can we do this in this routine.

Really appreciate!

Dear Lingling,

The UserTwrLd() routine of FAST v7 requires that you return the load (per unit length) at the tower node that is passed to the routine (separately for each node). The global displacement (X,Y,Z) and rotation (about X,Y,Z, assuming small angles) of the node are passed to the routine.

The example UserTwrLd() routine available in the UserTwrLd_DS.f90 file I attached in this forum topic above, interpolates a depth-varying spring to the actual depth of the node. You could apply the same principle for a depth-varying p-y curve.

I hope that helps.

Best regards,

Dear Jason,

Thank you for response!

I am trying to understand the way FAST pass the load and displacement of tower node to subroutine UserTwrLd.

If now FAST is passing load and displacement of JNode - TwrFt, X and I change the original format of source you give me:

zi ( 1) = -88.3 zi ( 2) = -88 zi ( 3) = -87 ..... Spring ( 1) = 1978166000 Spring ( 2) = 1944235000 Spring ( 3) = 1910304000 .....

zi ( 1) = -88.3 zi ( 2) = -88 zi ( 3) = -87 ..... P( 1) = a*TANH(b*X(1)) P( 2) = c*TANH(d*X(1)) P( 3) = e*TANH(f*X(1)) ...

, where a,b,c,d,e,and f are some constants for sand p-y curve relation and P is the soil resistance.

If JNode which FAST is current passing is located at zi ( 3) = -87 and, X(1) in my source equation is becoming the 1-direction displacement of JNode of tower (located at zi ( 3) = -87). And as a matter of fact, all the soil resistance P(1), P(2), P(3)… at different location of tower will be related to the 1-direction displacement of JNode of tower (located at zi ( 3) = -87).

Even though all the soil resistance P(1), P(2), P(3)… at different location of tower will be related to the 1-direction displacement of JNode of tower (located at zi ( 3) = -87), FAST is only using this p(3)-X(1) relation for this zi (3) Node calculation which is physically correct. By using InterpStp we can interpolate this soil resistance (right now is located at zi (3)=-87) as you did in your source Stff(J) and use Stff(JNode) to only obtain TwrFt at zi (3)=-87.

So if some other JNodes information are passing by, the p-y curves in this source will all be related to the displacement of this JNode but only the force TwrLd at JNode is determined by using the only correct p-y relation and that is the way to recompile it.

Hope I am doing things correctly.

Really appreciate!

Dear Lingling,

I’m not following everything you’ve said, but I would think it would be preferable to evaluate the force only at the given node. To do this, I suggest you interpolate the p-y curves to get the specific p-y curve at each node (depth) of the tower. For example, at model initialization you’d calculate p-y constants, a(J) and b(J), for each node (depth) of the tower (J = 1 to TwrNodes). Then, when UserPtfmLd() is called for JNode at every time step, you could calculate the load in the 1-direction at JNode as follows:

TwrFt(1) = -a(JNode)*TANH( b(JNode)*X(1) )

I hope that helps.

Best regards,

Dear Jason,

I am sorry for the confusion.

Could you please help to check the attached two documents when you get a chance:

  1. Please change source p-y.txt to source p-y.xlsx file. Inside this Excel document, the first sheet (dense sand (p-y)) is to pre-calculate soil resistance p = atanh(bX(1) at various depths. The second sheet (FAST (p-y)) is to organize pre-calculated p-y relation into FAST format.

  2. HydroCalc.txt is the FAST source I have compiled in order to incorporate p-y curves. I have made some changes at:

Line 5128 (Comment)
Line 5153
Line 5165 - 5287
Line 5289 (Comment)
Line 5343 (I dropped X(K), since I have multiplied this term in my p-y relation)

Really appreciate!
HydroCalc.txt (298 KB)
source p-y.txt (31.1 KB)

Dear Lingling,

I’m not sure what you want me to check in “source p-y.xlsx”. I can’t really confirm how you derived the p-y curves.

I don’t agree with how you’ve implemented the p-y curve within routine UserTwrLd(). As I explained in my prior post, I think you’d want line 5343 to be something like TwrFt(K) = TwrFt(K) - a(JNode)*TANH( b(JNode)*X(K) ), where a() and b() are defined at initialilization (within the FirstPas IF-THEN test).

Best regards,

Dear Mr Jonkman,

In my Master Thesis I am investigating the different parameters that can induce variability in the dynamic behavior of an onshore 5MW NREL Reference WT. I would like to ask some questions regarding Soil Structure Interaction for the case of an onshore WT. Mainly, I am working with FAST version 8, although I have also investigated how version 7 works as well as the services that it can offer and are not yet integrated in version 8 (i.e. linearization). I have read plenty of the past topics, but still have some points unclear.

As already mentioned in various topics, in order to model boundary conditions other than a cantilevered beam, it is possible to calculate the correct mode shapes of these boundary conditions and run the analyses. First allow me to sum up the different approaches, as I have organised them so far based on my findings:

1) Apparent Fixity (AF)
Applicable in FAST v.8.00. Within SubDyn, increase the length of the substructure (i.e. a Monopile), using a fictitious length and respective effective cross-sectional properties in the additional elements that have been introduced. Execute the stand-alone version of SubDyn and from the summary file extract the Mass and Stiffness matrices of the substructure. Plug them in BModes in the position of hydro_M and hydro_K. Change (?) the file that contains the distributed tower properties (due to the fictitious elements that belong to the AF length). Execute BModes, obtain from ModeShapePolyfitting.xls the correct polynomial coefficients of the mode shapes and run the analysis in FAST.

2) Coupled Springs (CS)
Only applicable in FAST v.7.00.
Follow similar procedure to generate altered mode shapes. Introduce influence of soil-structure interaction through a user-defined routine (file UserPtfmLd.f90) which allows implementing a discrete load at the platform reference point. 6x6 Stiffness and Damping matrices have to be defined in the routine. Then, recompilation of FAST is required.

3) Distributed Springs (DS)
Only applicable in FAST v.7.00
Follow similar procedure to generate altered mode shapes. Introduce influence of soil-structure interaction through user-defined routine (UserTwrLd.f90) which allows to implement distributed loads (per unit length) along the flexible tower. Might be possible to incorporate p-y curves. UserTwrLd interpolates a depth varying spring to the actual depth of the node. In the same sense the p-y curves can be interpolated in order to get a separate curve for its node. Also recompilation of FAST is required.


  1. When generating altered mode shapes in BModes, what is the role of the platform mass and (mass_pform) the platform mass inertia 3X3 matrix (i_matrix_pform)? I assume that in an onshore WT there is no platform; however, SubDyn does generate the concentrated total mass of the substructure and the coordinates of its center of mass. In another post you have mentioned that the platform mass inertia matrix as well hydro_M and hydro_K can be set as zero matrices for a monopile. How is that possible? When the tower boundary conditions in BModes are set to free-free (2-------hub_conn), the derived frequencies depend strongly on hydro_M and hydro_K matrices.

  2. Regarding the AF approach, should the AF length be added to the “draft” parameter in BModes?

  3. Since there is no mooring system in an onshore WT, I suppose that the Mooring-system 6X6 stiffness matrix (mooring_K) has to be set to zero. This influences quite significantly the frequencies generated by BModes (they decrease). Is this procedure correct or should this stiffness reduction be somehow made up for?

  4. I noticed that the Stiffness matrix within UserPtfmLd.f90 is the same with mooring_K in BModes and that the distributed springs in UserTwrLd.f90 are the same with the “Distributed elastic stiffness per unit length along a flexible portion of the tower length” in BModes. Could it be correct, that, by only changing the mode shapes when implementing the CS or DS models, we can substitute the respective subroutines and avoid the whole recompilation process? Isn’t it wrong on the contrary to do both and duplicate the effect through BModes and within the subroutines?

  5. Where exactly are the points in “Distributed elastic stiffness per unit length along a flexible portion of the tower length” defined? The user can set the number of the desired points and the respective spring stiffness, but how is their exact location along the monopile defined?

I apologize for the long text! I greatly appreciate your useful advice!

Best Regards,

Athanasios Kontis
MSc Civil Engineering Student, ETH Zurich

Dear Athanasios,

Your understanding of the AF, CS, and DS models is correct, except that in AF for FAST v8 with SubDyn, the file of distributed tower properties in BModes should only contain the “real tower”, and so, should be unaffected by AF; instead, the AF is accounted for in the hydro_M and hydro_K matrices taken from the SubDyn summary file and specified in the BModes primary input file.

To answer your direct questions:

  1. When generating mode shapes from BModes for use in FAST v8, cm_pform, mass_pform, and i_matrix_pform from BModes should match the corresponding platform mass/CM/inertia from the ElastoDyn primary input file (PtfmCMxt, PtfmCMyt, PtfmCMzt, PtfmMass, PtfmRIner, PtfmPIner, PtfmYIner). These can all be set to zero if that is appropriate for the model i.e. when there is no lumped mass/inertia between the tower and substructure. In BModes, hydro_M and hydro_K should not be zero unless a fixed tower-base boundary condition is enforced.

  2. No. Again, when using BModes to generate the tower mode shapes for modeling the AF approach with FAST v8, the tower in BModes should be the “real tower”, the AF (fictitious) section should be modeled in SubDyn and its effect should be represented through the hydro_M and hydro_K matrices in BModes. Input ref_msl from BModes should match PtfmRefzt from FAST and input draft from BModes functions similar to TowerBsHt in FAST (although they should have opposite signs).

  3. I’m not sure I understand your question, but hydro_K and mooring_K are simply summed in BModes, so as long as there sum is correct, BModes should function as expected.

  4. I’m not sure I understand your question.

  5. I’ve never used this feature of BModes, and so, cannot comment. Perhaps another of user of BModes can respond.

Best regards,