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 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
.....
to
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.
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:
Could you please help to check the attached two documents when you get a chance:
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.
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)
…
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).
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.
Questions:
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.
Regarding the AF approach, should the AF length be added to the “draft” parameter in BModes?
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?
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?
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
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:
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.
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).
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.
I’m not sure I understand your question.
I’ve never used this feature of BModes, and so, cannot comment. Perhaps another of user of BModes can respond.
Thank you for your input. I understand your point concerning the AF approach.
I will try to clarify my 4th question and through that also explain my strategy on implementing a Distributed Springs model. I am a bit unsure of the procedure I follow and I would really appreciate your opinion.
Generally, I try to avoid implementing springs through UserTwrLd.f90 and recompiling FAST, as in that case I would have to work with version 7. My experiments are completely based on the example presented in “OC3-Derivation and Description of the Soil-Pile-Interaction Models”; meaning I use the same WT, soil profile, properties etc, with the only difference being that I examine an onshore WT. The steps that I follow while working with version 8 are the following:
The penetration length of the monopile in the above example is 36m. I model the monopile in SubDyn from 0.0 m to -36.0 using appropriate number of elements and I run it in standalone mode to obtain the hydro_M and hydro_K matrices.
Generate new mode shapes. Within BModes I change the distributed tower data, the boundary conditions of tower to “free-free” (hub_conn=2), I plug in hydro_M and hydro_K matrices and set the mooring_K matrix equal to zero. At this point I also define the stiffness of the linear springs along the pile via “n_secs_k_distr: 36”, while the values of these stiffnesses are taken from an excel sheet that was referring to the experiments on the aforementioned paper and was uploaded in another topic in this forum. So my question was referring to this exact point: Does the implementation of these stiffnesses in BModes (and calculating the respective mode shapes) substitute the procedure of defining springs in UserTwrLd.f90 and recompiling FAST version 7?
Another important question refers again to the “draft” parameter. In my experiments I have set “draft=36m”. From your previous answer I suppose this is wrong since modeling the substructure in SubDyn and including it in the analysis by FAST accounts for this effect. However, what confuses me and makes me doubt of this statement, is that if I set draft=0 and execute BModes, the natural frequencies that I receive are higher than those of the fixed-base WT (i.e 1SS 0.57 Hz, 1FA 0.58 Hz), because of the implementation of 36 additional springs.
Then I generate the mode shape coefficients in ModeShapePolyfitting.xls and plug them in the tower data.
I activate 4 out of 6 Platform DOFs (heave and yaw deactivated), I activate SubDyn and run FAST.
Could you please share your opinion on the process I described above? I would really appreciate it.
Generally I post process my results through the PSDs of the tip acceleration of blade and tower. Comparing my experiments between a fixed base WT, an AF approach and a DS model I obtained the following frequency-peak results for the tower:
Fixed base:
1st peak: 0.315 Hz 2nd peak: 2.94 Hz
AF Approach:
1st peak: 0.273 Hz 2nd peak: 2.46 Hz
DS Model
1st peak: 0.273 Hz 2nd peak:1.64 Hz
I guess the reduction of the 1st frequency is reasonable and it’s encouraging that there is such a good agreement between AF and DS models. Where is the 2nd frequency expected though? Is it expected that the frequency peaks of the 2 models are different or does it mean that one of them is likely to have a problem?
In BModes, presumably you’ve only modeled the tower (above the ground) as flexible because you’re considering the flexibility of the pile (below ground) to be represented by hydro_M and hydro_K matrices. Thus, input draft in BModes should be set to zero. I’m not fully aware of how BModes functions, but I’d be surprised if you could apply linear springs to the section of the beam below ground when that section is not actually represented as a beam in BModes.
Adding springs to BModes to influence the tower mode shapes used in FAST is not a substitute for adding springs directly to the FAST model. You’d also have to modify FAST to directly model the springs.
The bottom of the pile should not be fixed (as it would be in SubDyn) when implementing a distributed-springs model.
So essentially what I understand is that the use of SubDyn in this case is out of the question, since the most recent version has to have at least one constrained node at the bottom. How does the following concept for the DS model sound to you:
calculate the flexible tower BModes, by setting hub_conn=2, leaving hydro_K, hydro_M and mooring_K matrices to zero and setting the tower draft=36m. Essentially leaving the substructure subprogram behind.
since adding springs to BModes does not substitute the direct implementation of springs in FAST, additionally determine this springs in the routine UserTwrLd.f90 and recompile FAST (which requires working with version 7.00)
Yes, that sounds like the correct approach for FAST v7. Just to be clear, when the draft is set to 36 m in BModes, the distributed tower properties should also include the properties of the pile below ground.
Your previous posts in this topic have been really helpfull, thank you, but there is still a few things I find unclear about updating the mode shapes.
I am modelling the 5MW OWT in FAST v8, and used the AF approach in SubDyn. I used test19.fst as a base. What I have done is:
In the SubDyn file:
The calculated AF length is 17.52 m, so I added another joint, located -37.52 m below MSL, and assigned the required properties.
I removed the contstraints of the mudline joint and constrained the new joint in all 6 DOFs.
In BModes:
In the BModes input file CS_monopile.bmi (provided by NREL), which is used for modelling the 5MW OWT monopile for rigid foundation, hydro_M and hydro_K are zero, while mooring_K contains values. To obtain the new mode shapes for the AF model, I replaced hydro_K and hydro_M with Kbbt and Mbbt matrices from the SubDyn summary file, while setting all mooring_K values to zero. This resulted in natural frequencies from BModes reasonably lower than for the rigid model (Mode 1 frequency = 0.189187 Hz for AF model vs 0.242302 Hz for rigid model).
Then I used ModeShapePolyfitting.xls to get the correct mode shape inputs to ElastoDyn, before running FAST for the AF model.
My question is: Is this the correct approach? I didn’t change anything else or any other input files for the test19.fst.
Something that confuses me is that Kbbt and Mbbt from the SubDyn summary file for the rigid model (without AF method) is not zero and not the same as mooring_K, while hydro_K and hydro_M is zero. If this is the correct approach, shouldn’t Kbbt and Mbbt from SubDyn for the rigid model be the same as hydro_K+moorinig_K and hydro_M in the Bmodes input file?
The CS_Monopile.bmi file that NREL has shared was developed for the old FAST v7 model of the OC3-monopile with flexible foundation represented via the coupled springs (CS) model. When employing the CS model for a monopile in FAST v7, the entire tower + monopile (down to the seabed) was modeled as the “tower” (SubDyn did not exist as a separate module), so, the whole tower + monopile was modeled in BModes as well, and the 6x6 stiffness matrix for the CS model was specified as mooring_K. For a rigid foundation, no matrices would need to be specified in BModes and a clamped boundary condition would be set.
In FAST v8, the tower is modeled in the ElastoDyn module and the monopile is modeled in the SubDyn module. BModes is used to derive the mode shapes of the tower for the ElastoDyn module, so in BModes, the monopile (whether with a rigid foundation or a flexible foundation via the AF model) is represented simply as stiffness (hydro_K or mooring_K) and mass (hydro_M) matrices.
Thank you very much for your response, I really appreciate it.
I see now that I misunderstood the CS_monopile files, of course they are for the coupled springs model.
Can I please ask you a few more questions just to be sure I have got eveything right?
I started with the CS_monopile.bmi input file and made the following changes for the AF model:
Changed Mooring_K to 0, and replaced hydro_K and hydro_M with the matrices from SubDyn, as explained in my last post.
Changed draft from 20 to 10 and ref_msl from 20 to 10.
Changed the ‘CS_Monopile_tower_secs.dat’ file to only contain values for the tower and not the substructure, i.e. I removed the two first lines of the original file (which makes the remaining properties the same as in the NRELOffshrBsline5MW_OC3Monopile_ElastoDyn_Tower.dat file). In addition, I changed flp_iner and edge_iner to very small numbers and tor_stff and axial_stff to very high numbers to avoid axial, torsion and shear modes.
Then I executed BModes and got the following results for the first modes (for the AF model):
Mode 1: 0.185392 Hz, 1st fore-aft
Mode 2: 0.185914 Hz, 1st side-to-side
Mode 3: 1.099990 Hz, 2nd side-to-side
Mode 4: 1.19299 Hz, 2nd fore-aft
Just to check if my approach was correct I tried to do the same for the rigid foundation model, run test19.fst and copy Mbbt and Kbbt from the SubDyn summary file into Hydro_M and Hydro_K, put Mooring_K to zero and ran BModes with the same changes as above. This resulted in the following results:
Mode 1: 0.205967, 1st side-to-side
Mode 2: 0.207538, 1st fore-aft
Mode 3: 1.289402, 2nd side-to-side
Mode 4: 1.450130 2nd fore-aft
I plotted the mode shapes and compared them with the mode shapes in the ElastoDyn input file provided in the FAST archieve for the ridig foundation model (NRELOffshrBsline5MW_OC3Monopile_ElastoDyn_Tower), and noticed that they differed slightly.
Is this the approach used to get the mode shapes input to ElastoDyn? Or how was it done? Is there anywhere I can find the BModes input/output file for the NREL 5MW OWT tower for use in Fast v8?
Does my natural frequencies for the tower with and withouth AF method seem correct?
My first thought was to compare the natural frequencies with the frequencies from the CS_monopile Bmodes output file (mode 1 = 0.242302 Hz), but these are for the whole tower+substructure, while the natural frequencies in my approach are for only the tower. Are these frequencies comparable?
Why is draft and ref_msl 10 m and not -10? In the BModes file used with FAST v7 the default value is 20m, because the tower base is 20 m below MSL. In FAST v8, where you’re only supposed to calculate mode shapes for the tower, starting from 10m above MSL, shouldn’t draft then be set to minus 10 meters?
Is there anything else you think I have misunderstood?
You’re overall approach sounds correct, but your question (4) is spot-on. Instead of setting draft = ref_msl = 10, you should be setting draft = ref_msl = -10 (these inputs are defined positive downward in BModes).
Here are my answers to your direct questions:
Actually, I don’t have the BModes input files used to derive the tower mode shapes for Test19. These were developed by a post-doc, who is no longer employed at NREL. Perhaps changing draft and ref_msl as indicated above will give you closer results.
Your natural frequencies look low, but I suspect this is because your tower is 20-m longer that it should be (based on the incorrect draft and ref_msl).
Actually, the frequencies predicted by BModes will be for the tower + substructure in your case, its just that the tower is modeled with finite elements while the substructure is modeled as mass and stiffness matrices, so, the shapes predicted by BModes will only be for the tower (with a moveable base).
See above.
Your specification of draft and ref_msl is the only issue I can see.
The frequencies are higher, but shouldn’t the fore-aft and side-side modes frequencies be more similar? I would expect Mode 1 (1st fore-aft) and Mode 2 (1st side-side) to have almost identical frequencies, and the same with Mode 3 and 4.
I can’t find an explanation to this, can you please tell me what you think?
I have attached my BModes input files for the rigid and the AF model along with the SubDyn input file for the AF model (after new calculations of soil stiffness I changed the AF length to 14.2784m).
Also, I have seen in that users of this forum ask for “the unofficial FAST theory manual” and user’s guide for FAST v8, could you please send me a copy as well? My email is: ingbylok@gmail.com
A new question has come up that I hope you can help me with.
I am wondering why the natural frequencies from BModes for the original model (rigid foundation, no modification done in any of the FAST input files provided by NREL) is not expected to be the same as the natural frequencies given in “Definition of a 5-MW Reference Wind Turbine for Offshore System Development”. What is the difference between these “full system natural frequencies” and the ones from BModes? The first fore-aft frequency calculated from BModes is 0.2804 Hz, while the one given in the report is 0.3240 Hz.
I computed the PSD from the fore-aft tower top displacement (excited the structure with white noise wave spectrum only, no wind, as explained in this forum post: AeroDyn Feedback), the PSD figure is attached, the rigid (original) foundation model vs the apparent fixity model. The frequency peaks seem to be approximately similar to the fore-aft natural frequencies from BModes, but not exactly. Does this PSD seem correct? What have I misunderstood about the “full system natural frequencies?”, what do they represent and why are they not the ones present in the PSD plots?"
In general, BModes and FAST may predict different tower natural frequencies because the rotor-nacelle assembly (RNA) is modeled rigidly in BModes, but flexible in FAST (the coupling between the rotor and tower is more important for the second tower-bending modes than the first). If you disable the degrees of freedom in the rotor, drivetrain, and nacelle of the FAST model, and set up the BModes and FAST models consistently, I would expect that you’d get quite similar predictions of tower natural frequencies between the tools.
In your case, though, the first tower-bending natural frequency of 0.32 Hz reported in the NREL 5-MW specifications report is for the land-based wind turbine. If you are exciting your model with white noise waves, I suspect you’re using the NREL 5-MW turbine atop the OC3-monopile, which has a first support-structure natural frequency of 0.28 Hz, which matches well what you are reporting from BModes and your PSD.
Yes, I am using the NREL 5MW offshore monopile turbine, and had mistakenly though that the natural frequencies in the specifications report were for the offshore turbine, so this clarifies a lot. I am satisfied with my results now!
