Modeling of linear foundation for onshore turbine

Hi everybody,
I am trying to convert a model from Bladed to OpenFAST as accurate as possible. In Bladed the tower base can be either coupled to the inertial frame or to a specified mass/inertia via a single linear spring per DOF with specified stiffness. I would like to reproduce the later option with OpenFAST.

In this post ([url]Foundation Stiffness and damping in FAST v8]) it is mentioned that this could be done using HydroDyn without the “hydro”-parts. Is that still the recommended way to go? If yes, could some one give me an example file? I am not sure how to set some of the parameters:
[] Do WtrDens, WtrDpth and MSL2SWL have an effect, if WaveMod, CurrMod and PotMod are zero?
[
] Should the platform force flags be true or false?
[] Can I set NAxCoef, NJoints, NPropSets, NCoefDpth, NMembers, NFillGroups to zero?
[
] What should I choose for the simple hydro dynamic coefficients (model 1)

Thanks already for you help!

Dear Jens,

I gather from your question that you are modeling a land-based wind turbine; is that correct? If so, the simplest solution may be to enable the user-specified external platform (ExtPtfm) by setting CompSub = 2 and CompHydro = 0, which is briefly documented here: openfast.readthedocs.io/en/mast … nfast.html.

Best regards,

Hi Jason,
thanks for the quick reply.
Indeed I am modeling an onshore turbine. I will edit the title accordingly.
Is the file format for the SubFile the same as when CompSub=1? Or is there any place I can find an example file (already looked in the matlab-toolbox and r-test repos)?

BR

Dear Jens,

Good question; no, the file format is specific to ExtPtfm. I looked a bit and can’t seem to find a sample input file. But looking at ExtPtfm_MCKF.f90/ReadPrimaryFile(), the format of the input file looks quite clear:

Lines Description 1 File Header: External Platform MCKF Matrices 2 Section Header: Mass Matrix 3-8 PtfmAM 9 Section Header: Damping Matrix 10-15 Damp 16 Section Header: Stiffness Matrix 17-22 Stff 23 Section Header: Loads time-history 24 Loads time-history table channel names, i.e. Time PtfmFt_FX PtfmFT_FY PtfmFT_FZ PtfmFt_MX PtfmFT_MY PtfmFT_MZ 25 Loads time-history table channel units, i.e. (s) (N) (N) (N) (Nm) (Nm) (Nm) 26-end PtfmFt_t PtfmFt
I hope that helps.

Best regards,

Great, thanks I’ll try this and report back how it works.

Here is my promised report. From what I saw so far it seems to work as expected.
In the main fst file I set:

          2   CompSub         - Compute sub-structural dynamics (switch) {0=None; 1=SubDyn; 2=External Platform 

and

"sub_file_path"      SubFile         - Name of file containing sub-structural input parameters (quoted string)

In the ED file I set, e.g.

True         PtfmSgDOF   - Platform horizontal surge translation DOF (flag)
True         PtfmSwDOF   - Platform horizontal sway translation DOF (flag)

and I have a SubExtPtfm file like this:

------- External Platform MCKF Matrices INPUT FILE -----------------------------
------- MASS MATRIX ------------------------------------------------------------
    1.000000E+00     0.000000E+00     0.000000E+00     0.000000E+00     0.000000E+00     0.000000E+00 
    0.000000E+00     1.000000E+00     0.000000E+00     0.000000E+00     0.000000E+00     0.000000E+00 
    0.000000E+00     0.000000E+00     1.000000E+00     0.000000E+00     0.000000E+00     0.000000E+00 
    0.000000E+00     0.000000E+00     0.000000E+00     1.300000E+07     0.000000E+00     0.000000E+00 
    0.000000E+00     0.000000E+00     0.000000E+00     0.000000E+00     1.300000E+07     0.000000E+00 
    0.000000E+00     0.000000E+00     0.000000E+00     0.000000E+00     0.000000E+00     1.000000E+00 
------- DAMPING MATRIX ---------------------------------------------------------
    1.000000E-02     0.000000E+00     0.000000E+00     0.000000E+00     0.000000E+00     0.000000E+00 
    0.000000E+00     1.000000E-02     0.000000E+00     0.000000E+00     0.000000E+00     0.000000E+00 
    0.000000E+00     0.000000E+00     1.000000E-02     0.000000E+00     0.000000E+00     0.000000E+00 
    0.000000E+00     0.000000E+00     0.000000E+00     7.211103E+06     0.000000E+00     0.000000E+00 
    0.000000E+00     0.000000E+00     0.000000E+00     0.000000E+00     7.211103E+06     0.000000E+00 
    0.000000E+00     0.000000E+00     0.000000E+00     0.000000E+00     0.000000E+00     1.000000E-02 
------- STIFFNESS MATRIX -------------------------------------------------------
    1.000000E+00     0.000000E+00     0.000000E+00     0.000000E+00     0.000000E+00     0.000000E+00 
    0.000000E+00     1.000000E+00     0.000000E+00     0.000000E+00     0.000000E+00     0.000000E+00 
    0.000000E+00     0.000000E+00     1.000000E+00     0.000000E+00     0.000000E+00     0.000000E+00 
    0.000000E+00     0.000000E+00     0.000000E+00     4.000000E+10     0.000000E+00     0.000000E+00 
    0.000000E+00     0.000000E+00     0.000000E+00     0.000000E+00     4.000000E+10     0.000000E+00 
    0.000000E+00     0.000000E+00     0.000000E+00     0.000000E+00     0.000000E+00     1.000000E+00 
------- LOADS TIME HISTORY -----------------------------------------------------
Time PtfmFt_FX PtfmFT_FY PtfmFT_FZ PtfmFt_MX PtfmFT_MY PtfmFT_MZ
 (s)  (N)       (N)       (N)       (Nm)      (Nm)      (Nm)

This is written by a MATLAB function which essentially looks like this:

M= diag([m m 1 I1 I2 1]);
K= diag([K_trans K_trans 1 K_rot K_rot 1]); 
D= diag(0.005*2*sqrt(diag(K).*diag(M)));


fprintf(fid, '------- External Platform MCKF Matrices INPUT FILE -----------------------------\n');
fprintf(fid, '------- MASS MATRIX ------------------------------------------------------------\n');
for row= 1:6
    for col= 1:6
        fprintf(fid, '   % 12.6E ', M(row, col));
    end
    fprintf(fid, '\n');
end
fprintf(fid, '------- DAMPING MATRIX ---------------------------------------------------------\n');
for row= 1:6
    for col= 1:6
        fprintf(fid, '   % 12.6E ', D(row, col));
    end
    fprintf(fid, '\n');
end
fprintf(fid, '------- STIFFNESS MATRIX -------------------------------------------------------\n');
for row= 1:6
    for col= 1:6
        fprintf(fid, '   % 12.6E ', K(row, col));
    end
    fprintf(fid, '\n');
end
fprintf(fid, '------- LOADS TIME HISTORY -----------------------------------------------------\n');
fprintf(fid, 'Time PtfmFt_FX PtfmFT_FY PtfmFT_FZ PtfmFt_MX PtfmFT_MY PtfmFT_MZ\n');
fprintf(fid, ' (s)  (N)       (N)       (N)       (Nm)      (Nm)      (Nm)\n');

Thanks!

Hi Jens ,
i would like to ask you about the calculation inside the 2 Matrix(Damping, Stiffness).
could you please explain me more about this results .
Best regards ,

Hi meriem,
I set the matrices as follows:

  • all matrices are diagonal because I didn’t have more detailed data, but diagonal is usually a good choice in this case.
  • the damping matrix must be a real physical damping coefficient “d” (never mind the capital variable name in the script) but I only had a damping ratio “D” given that was 0.5% (0.0005). The calculation of damping coefficient from damping ratio is d=2Dsqrt(m*k). If you are not familiar with this relationship you can look it up in the wikipedia.

Finally if you don’t have a mass given, this means your foundation is coupled to the inertial frame. FAST does not allow this directly but you can achieve almost the same effect by setting the mass and moments of inertia to very high values. Ideally you should choose the mass and moment of inertia of the earth as values but that might cause numerical problems. Usually values about 10 to 50 times those of the total turbine should emulate an inertial foundation well enough.

Here is another question and answer:
Q:
but how you choose the damping ratio as 0,0005 in some doc i found that they took as 0,15 ?
why all the values of diag are 1 just 2 of them are different? (mass ans stiffness) …

A:
0.5% damping is quite usual. But 0.15% (or is it 15%?) may also work. You have to run some simulations to see what suites you best.
I only had values for translation and rotation about the x- and y-axes given, so I set the z-axes values to one. But these values only make sense when they are not being used, i.e. the corresponding degrees of freedom are set to false in the ElastoDyn config file. Otherwise such a low stiffness would make the turbine sink miles into the ground :frowning:

And another question I received by mail:
Q:
The mass moment of inertia I1,I2 How you choose their values ?

A:
All masses, moments of inertia and stiffnesses I worked with came from a Bladed model. If you don’t have such reference values, you might want to make the foundation like 10 times more inert that your turbine. For the moments of inertia you could take the mass of the rotor nacelle assembly and multiply it by the square of the tower hight as an approximation of the turbine values.

Hi Jens,

I am also converting a Bladed model to OpenFAST. I’m just wondering how you dealt with the missing DOFs in OpenFAST? That is: tower torsion, nacelle tilt, and main shaft bending? I found a workaround for tower torsion here (Tower Eigenfrequencies of NREL 5MW Turbine - #60 by Jason.Jonkman) but I’m wondering about the other two.

Thanks!

Best,
Drew Gertz
Principal Consultant
Northwind Engineering OÜ
www.northwindengineering.com

Hi Drew,
I just plain ignored these missing DOFs. In typical turbines they shouldn’t make much of a difference. The biggest influence is probably that of the tower torsion. I didn’t read the topic linked by you but I assume it suggests to emulate tower torsion with a yaw stiffness which should work quite well.
If a very good comparison must be achieved, those DOFs must be disabled in Bladed.

Hi,
My apologies for the late reply, thank you for your detailed and helpful explanation.
i’am studying onshore wind turbine and i want to know the effect of the foundation on the tower loads , i would like to ask about how could i calculate the eigenfrequencies (with the consideration of the foundation effect) , and how can i set the rigidity parameter ,could i set this parameters ( PtfmSgDOF ,PtfmSwDOF,PtfmHvDOF) in Elastodyn True,and if could i there’s some conditions? or it’s not possible ?
Thank you.
Kind regards,

Dear Meriem,

Are you using OpenFAST with a the external platform (CompSub = 2) option? The linearization capability of OpenFAST does not currently support this option, although you could derive the eigenfrequencies by post-processing of time-domain responses, as has been discussed on this forum.

NREL has been working on modeling foundation flexibility and soil-structure interaction in OpenFAST by both including a 6x6 stiffness matrix (coupled springs model) at joints in SubDyn and by introducing a new module, SoilDyn. SoilDyn provides an option for a simple 6x6 stiffness matrix, but also includes options for interfacing to the REDWIN DLLs developed by NGI (to model soil-structure interaction via a superelement, including soil hysteresis), as well as and a placeholder for future implementation of a p-y and t-z curves. These are not currently included in a public release, but are being developed in branches on forks of OpenFAST. And these developments also support full-system linearization capability for direct eigenanalysis. More information will come soon.

Best regards,

Hi Jason,
thank you for your fast answer,
Yes, I’m using CompSub=2.
Sorry if I were not that much clear at first, but I would like to check the effect of the foundation on the tower loads by giving foundation’s mass, inertia and stiffness. To start a simulation with OpenFAST I need to give different shape modes of the tower based on the natural frequencies, which are calculated with the help of Bmodes v1.03.01 (Elastodyn_tower.dat). So far I simulated wind turbines with a rigid foundation, therefore I didn’t need to consider foundation effects on the tower natural frequency, which I need to do now! My first question is how can I consider this effect of the foundation on the tower’s natural frequency?
The other point is I would like to compare the calculated loads based on a rigid foundation (CompSub=0) with a modelled foundation (CompSub=2). I don’t need to calculate the loads on the foundation and as I know the platform degrees of freedom ( in Elastodyn.dat) only support the offshore foundations.
The question is, shall I turn on the degrees of freedom of the platform by setting them “True” or not? and does OpenFAST consider the given foundation parameters by the external file (External Platform MCKF) to calculate the tower’s loads when these degree of freedoms are “false” or not?

False PtfmSgDOF - Platform horizontal surge translation DOF (flag)
False PtfmSwDOF - Platform horizontal sway translation DOF (flag)
False PtfmHvDOF - Platform vertical heave translation DOF (flag)
False PtfmRDOF - Platform roll tilt rotation DOF (flag)
False PtfmPDOF - Platform pitch tilt rotation DOF (flag)
False PtfmYDOF - Platform yaw rotation DOF (flag)

Kind regards,
Meriem Mbarek

Dear Meriem,

Just a few comments:

  • To model a flexible foundation, you can use CompSub = 2 or wait until the new SoilDyn/SubDyn capabilities are available. Regardless, if you intend to use CompSub = 2 and model the tower in ElastoDyn, then you should specify the mass/stiffness/damping of the foundation in ExtPtfm and modify the tower mode shapes in ElastoDyn based on the mass/stiffness of the foundation (tower-base boundary condition); BModes_JJ can be used to calculate the tower mode shapes. Any change to the mass/stiffness of the foundation in ExtPtfm or to the tower mode shapes in ElastoDyn will impact the tower natural frequency.
  • CompSub = 2 and the ElastoDyn platform DOFs work for either onshore or offshore wind turbines.
  • If all of the ElastoDyn platform DOFs are disabled, then ExtPtfm will not have any effect; ExtPtfm will only have an effect on the solution if at least one ElastoDyn platform DOF is enabled.

Best regards,

Dear Jason,

Thank you for your help.

best regards,

Dear Jason ,
Sorry, i would like ask about BModes_jj:
i tried first time without considering the foundation(hub_conn=1), and i have as a result for the natural frequency:
eigenvalue( 1) = 0.493762D+11 mode 1 frequency = 0.370346
When i consider the foundation i have as a resulat(hub_conn=3) :
eigenvalue( 1) = 0.436656D+05 mode 1 frequency = 0.000348
i would expect to find the natural frequency few percent lower than this value with the consideration of the interaction between foundation and the soil with the consideration of the rotation about horizontal axes wind turbine.
[b]====================== BModes v3.00 Main Input File ==================
NREL 5MW Tower

--------- General parameters ---------------------------------------------------------------------
true Echo Echo input file contents to *.echo file if true.
2 beam_type 1: blade, 2: tower (-)
0. romg: rotor speed, automatically set to zero for tower modal analysis (rpm)

  1.    romg_mult:  rotor speed muliplicative factor (-)
    

87.6 radius: rotor tip radius measured along coned blade axis, OR tower height above ground level [onshore] or MSL offshore
0. hub_rad: hub radius measured along coned blade axis OR tower rigid-base height (m)
0. precone: built-in precone angle, automatically set to zero for a tower (deg)
0. bl_thp: blade pitch setting, automatically set to zero for a tower (deg)
3 hub_conn: hub-to-blade or tower-base boundary condition [1: cantilevered; 2: free-free; 3: only axial and torsion constraints] (-)
20 modepr: number of modes to be printed (-)
t TabDelim (true: tab-delimited output tables; false: space-delimited tables)
f mid_node_tw (true: output twist at mid-node of elements; false: no mid-node outputs)

--------- Blade-tip or tower-top mass properties --------------------------------------------
3.500003109E+005 tip_mass blade-tip or tower-top mass (kg)
-0.4137754432 cm_loc tip-mass c.m. offset from the tower axis measured along x-tower axis (m)
1.9669893542 cm_axial tip-mass c.m. offset tower tip measures axially along the z axis (m)
4.370E7 ixx_tip blade lag mass moment of inertia about the tip-section x reference axis (kg-m^2)
2.353E7 iyy_tip blade flap mass moment of inertia about the tip-section y reference axis (kg-m^2)
2.542E7 izz_tip torsion mass moment of inertia about the tip-section z reference axis (kg-m^2)
0. ixy_tip cross product of inertia about x and y reference axes(kg-m^2)
1.169E6 izx_tip cross product of inertia about z and x reference axes(kg-m^2)
0. iyz_tip cross product of inertia about y and z reference axes(kg-m^2)

--------- Distributed-property identifiers --------------------------------------------------------
1 id_mat: material_type [1: isotropic; non-isotropic composites option not yet available]
‘CS_monopile_tower_secs.dat’ : sec_props_file name of beam section properties file (-)

Property scaling factors…
1.0 sec_mass_mult: mass density multiplier (-)
1.0 flp_iner_mult: blade flap or tower f-a inertia multiplier (-)
1.0 lag_iner_mult: blade lag or tower s-s inertia multiplier (-)
1.0 flp_stff_mult: blade flap or tower f-a bending stiffness multiplier (-)
1.0 edge_stff_mult: blade lag or tower s-s bending stiffness multiplier (-)
1.0 tor_stff_mult: torsion stiffness multiplier (-)
1.0 axial_stff_mult: axial stiffness multiplier (-)
1.0 cg_offst_mult: cg offset multiplier (-)
1.0 sc_offst_mult: shear center multiplier (-)
1.0 tc_offst_mult: tension center multiplier (-)

--------- Finite element discretization --------------------------------------------------
61 nselt: no of blade or tower elements (-)
Distance of element boundary nodes from blade or flexible-tower root (normalized wrt blade or tower length), el_loc()
0 0.003481894 0.010445682 0.017409471 0.024373259 0.031337047 0.038300836 0.045264624 0.052228412 0.059192201 0.066155989 0.073119777 0.080083565 0.087047354 0.094011142 0.10097493 0.107938719 0.114902507 0.121866295 0.128830084 0.135793872 0.13990 0.149721448 0.156685237 0.163649025 0.170612813 0.177576602 0.18454039 0.191504178 0.198467967 0.205431755 0.212395543 0.219359331 0.22632312 0.233286908 0.240250696 0.247214485 0.250696379 0.320334262 0.37971 0.424791072 0.45961 0.486635 0.51366 0.54068 0.5677 0.594715 0.62173 0.64875 0.67577 0.70279 0.72981 0.75683 0.78385 0.81087 0.83789 0.864905 0.89192 0.91894 0.94596 0.97298 1.0

--------- Properties of tower support subsystem (read only if beam_type is 2) ------------
0 tow_support: : aditional tower support [0: no additional support; 1: floating-platform or monopile with or without tension wires] (-)
0.0 draft : depth of tower base from the ground or the MSL (mean sea level) (m)
0.0 cm_pform : distance of platform c.m. below the MSL (m)
0 mass_pform : platform mass (kg)
Platform mass inertia 3X3 matrix (i_matrix_pform):
0. 0. 0.
0. 0. 0.
0. 0. 0.
0.0 ref_msl : distance of platform reference point below the MSL (m)
Platform-reference-point-referred hydrodynamic 6X6 matrix (hydro_M):
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 5.650000E+7 0. 0.
0. 0. 0. 0. 5.650000E+07 0.
0. 0. 0. 0. 0. 0.
Platform-reference-point-referred hydrodynamic 6X6 stiffness matrix (hydro_K):
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 4.356000E+10 0. 0.
0. 0. 0. 0. 4.356000E+10 0.
0. 0. 0. 0. 0. 0.
Mooring-system 6X6 stiffness matrix (mooring_K):
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.

Distributed (hydrodynamic) added-mass per unit length along a flexible portion of the tower length:
0 n_secs_m_distr: number of sections at which added mass per unit length is specified (-)
0. 0. : z_distr_m [row array of size n_added_m_pts; section locations wrt the flexible tower base over which distributed mass is specified] (m)
0. 0. : distr_m [row array of size n_added_m_pts; added distributed masses per unit length] (kg/m)

Distributed elastic stiffness per unit length along a flexible portion of the tower length:
0 n_secs_k_distr: number of points at which distributed stiffness per unit length is specified (-)
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 : z_distr_k [row array of size n_added_m_pts; section locations wrt the flexible tower base over which distributed stiffness is specified] (m)
595318000.0 1165208000 1129400000 1095553000 1059931000 1024493000 989209000 953643000 918718000 883287000 847803000 812541000 777187000 741870000 706616000 671440000 636229000 600957000 565919000 530470000 495081000 459574000 385327000 305479000 280059000 254125000 227500000 200112000 171927000 143115000 114173000 80184000 52237000 35561000 20912000 9000000 1156000 : distr_k [row array of size n_added_m_pts; distributed stiffness per unit length] (N/m^2)

Tension wires data
0 n_attachments: no of wire-attachment locations on tower [0: no tension wires] (-)
3 3 n_wires: no of wires attached at each location (must be 3 or higher) (-)
6 9 node_attach: node numbers of attacments location (node number must be more than 1 and less than nselt+2) (-)
0.e0 0.e0 wire_stfness: wire spring constant in each set (see users’ manual) (N/m)
0. 0. th_wire: angle of tension wires (wrt the horizontal ground plane) at each attachment point (deg)

END of Main Input File Data *********************************************************************
*************************************************************************************************[/b]

Best regards

Dear Meriem,

With hub_conn = 3, the tower-base (platform) point has 4 DOFs–surge, sway, roll, and pitch. It looks like you’ve specified the lumped inertia and stiffness in roll/pitch, but you have not specified any stiffness in surge and sway. This will result in two rigid-body modes (in surge and sway), each with zero frequency (or near-zero frequency due to numerical round off in the solution).

Best regards,