Onshore wind turbine: modeling of nonlinear foundation

Dear Jason,

If I do the evaluation at the platform level PtfmRDyi I get the following results:
1
If I do the evaluation at TwHt1RDyt I get the following results:
2
In both cases they are not right.

By fixing all the DOF’s at the platform I get the following results at TwHt1RDyt, which are in the expected range.
3

Thank you.

Best regards,
Mário

Dear @Mario.Sobrinho,

Regarding PtfmRDyi, this is the pitch rotation at the transition piece as output by ElastoDyn. From what I can tell, the turbine is simply falling over, suggesting SoilDyn is not providing a reasonable reaction moment with DLL_model = 1. The platform-pitch DOF calculation in ElastoDyn makes use of a moderate angle correction that loses considerable accuracy beyond about 20 degrees. So, I would not expect an accurate result after about 8 s into your simulation.

Best regards,

Dear @Mario.Sobrinho,

As commented before, I would run the simulation with the next settings:

  1. ElastoDyn: PtfmSgDOF = PtfmSwDOF = PtfmHvDOF = PtfmRDOF = PtfmPDOF = PtfmYDOF = True
  2. SubDyn: JointType = 1 and all DOFs constrained in the reaction joint except the rotation around Y axis. This is: RctTDXss = RctTDYss = RctTDZss = RctRDXss = RctRDZss = 1. RctRDYss = 0.
  3. OpenFAST input file: NumCrctn = 1 (one correction iteration).

Now, regarding SoilDyn it seems that is quite challenging with REDWIN model 1. So, as a first attempt to understand what is going on, I would propose you to use CalcOption 1 and define a linear stiffness (K55 position in the stiffness matrix). Hopefully, this will work as expected. Fingers crossed…

If you can successfully run with CalcOption 1, then you can switch to CalcOption 3 and try the REDWIN model 1 with a simple force/deflection relationship similar to what you used in CalcOption 1 (K55 stiffness). So we are adding complexity into the model step by step.

Roger

1 Like

Dear Jason, dear Roger,

Thank you for your replies.

Dear @Roger.Bergua, I ran a simulation considering the 3 points that you mentioned. In SoilDyn I used CalcOption 1 with a linear stiffness value in K55 as you suggested. Unfortunately I got the same results as before, i. e., model stops at the beginning.

I decided to lock all the DOF’s of the platform, except the PtfmPDOF. I did it for the ElastoDyn input file, as well as for the SubDyn input file. With this approach, I was able to run the model without any error message.
When I switch to CalcOption 3, the system fails as before. For this case, I used a linear relation moment/rotation with the same stiffness used in CalcOption 1.

I have one additional question related with the rotation output. By looking at the rotations at the platform (PtfmRDyi) and at the base of the tower (TwHt1RDyt), the numbers are very different, with the rotation at the platform strangely high, as you can see from the plots below. Even if the rotation of the platform was correct, I would expect a higher rotation in the tower due to some deflection and not the opposite. Do you have any idea why is this happening?
4

Thank you.

Best regards,
Mário

Dear @Mario.Sobrinho,

I’m surprised the model is not running with CalcOption = 1 in SoilDyn. Did you set the (5,5) element of the stiffness matrix in SoilDyn appropriately?

I would not recommend disabling all ElastoDyn platform DOFs except PtfmPDOF. This would eliminate all Guyan modes except for the interface pitch mode. Again, you should keep all of the ElastoDyn platform DOFs enabled and only lock the pitch DOF of the SubDyn reaction joint.

ElastoDyn output PtfmRDyi is the rotation of the platform reference point in the inertial (i) coordinate system; this angle could get large if the tower is unconstrained (as shown by your results). ElastoDyn output TwHt1RDyt is the rotation of the tower output node 1 relative to the platform in the tower-base coordinate system; I would expect this angle to remain small because the tower deflections should be small.

Best regards,

Dear Jason,

Yes, I believe I did it right. Please, see below:

------- SoilDyn v1.0 INPUT FILE ------------------------------------------------------------------------------
Example SoilDyn input file for stiffness matrix in the seabed.
--------------------------------------------------------------------------------------------------------------
      False                Echo           - Echo input data to <RootName>.ech (flag)
      "default"            DT             - Communication interval (s) (or "default")
      1                    CalcOption     - Calculation option  {1: Stiffness / Damping matrices, 2: P-Y curves [unavailable], 3: coupled REDWIN DLL}
================== Parameters for Stiffness / Damping matrices [used only for CalcOption = 1] ================
   0 0 -0.2                   Location (X,Y,Z)     - the location of interface point (only one point allowed)
Stiffness matrix (6x6)
   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00
   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00
   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00
   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00
   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00   0076013808   0.0000e+00
   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00
Damping ratio matrix (6x6)
   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00
   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00
   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00
   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00
   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00
   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00

If I understood correctly your answer related to the rotations, the TwHt1RDyt is the rotation relative to the rotated platform? The accumulated rotation at node 1 would be the PtfmRDyi + TwHt1RDyt. Is that correct?

Thank you.

Best regards,
Mário

Dear @Mario.Sobrinho,

Does the (5,5) element of the stiffness matrix actually have leading zeros? I’m not sure how this would be interpreted.

I agree with your calculation of the total rotation of ElastoDyn tower output node 1.

Best regards,

Dear Jason,

Same result :frowning:

Best regards,
Mário

Dear @Mario.Sobrinho,

How big (or small) is your wind turbine? The bending stiffness that you are defining is 7.6E7 Nm/rad. This value is certainly small. For example, the rated thrust force of the IEA-10.0-198-RWT is close to 1.5E6 N. Assuming a distance between the tower top and the seabed of 100 m, the bending moment at the seabed would be 1.5E8 Nm. This bending moment would make your monopile rotate 1.97 rad…

In the section 3.3.1 of the next technical report (https://www.nrel.gov/docs/fy21osti/79938.pdf), you can see the stiffness matrix used during the OC6 phase II project. The bending stiffness is 8.1E11 Nm/rad. This is 4 orders of magnitude stiffer than the value that you are using.

Finally, it is important to note that if you are using CalcOption 1 in SoilDyn you are just defining a stiffness matrix as the boundary condition for the SubDyn reaction joint. You can include that stiffness matrix in SoilDyn (as you are doing) or you can include the stiffness matrix in SubDyn directly. Maybe you want to start using SubDyn-only to ensure that your system is stable and add complexity later on. For reference, during the OC6 phase II project I used both approaches (SubDyn-only and SubDyn+SoilDyn CalcOption 1) and both approaches returned the exact same behavior (as expected). Below you can see the definition used for OC6 phase II when using the SubDyn-only approach. This stiffness matrix corresponds to the one shown in the section 3.3.1 that I commented before.

I hope that helps.

Roger

1 Like

Dear Roger,

Thank you for the shared information. I already tested for different range of stiffnesses but the result is always the same, even using the values you shared.
The WTG is with 100 m rotor diameter and 95 m hub high. I’m testing with a simple wind flow by using steady wind type with 8 m/s of wind speed (InflowWind).

Once the model starts I receive always this message for the nodal output. I don’t think this is related with the problem, but in any case I just would like to mention this to you.

Could be something related with the WTG configuration defined in the ElastoDyn?
I tested different coordinates for the note 2 in SubDyn and SoilDyn, but without any changes on the results. Below you can find a sketch of the WTG and foundation.

Thank you.

Best regards,
Mário

Dear @Mario.Sobrinho,

It’s quite surprising to me… I’m not sure if the issue that you are having is related to SubDyn/SoilDyn or with your numerical model somewhere else. Btw, the message about nodal outputs is normal. It’s just informing you that you didn’t provide information for the nodal output section (i.e., advanced sensors that only some users request).

I’m wondering if your model is working as expected when you use SubDyn with a clamp condition at the base reaction joint?

If I would be you, I would increase the complexity in the model one step at a time. I would not move to the next step (increase complexity) if things are not looking as expected. These are the steps that I would follow:

  1. Run the model with ElastoDyn (no SubDyn and no SoilDyn). Is the system working as expected? Does the bending moment at the bottom look good?
  2. Run the model including SubDyn with a clamp at the base. You should get equivalent results to point (1).
  3. Run the model with SubDyn. Remove the clamp and include a stiffness matrix in the diagonal terms that is quite stiff (not too stiff because otherwise the system will become ill-conditioned). With these settings you should reproduce the results from points (1) and (2).
  4. Run the model with SubDyn. Define the bending stiffness desired for the Y direction. Is the system rotating as expected? With the bending moment obtained in points (1-2-3) and the bending stiffness, you can anticipate how much the system should rotate at the base (approx. value)
  5. Run the model with SubDyn and SoilDyn. Remove the bending stiffness in SubDyn and include it in SoilDyn CAlcOption 1. This should return the exact same results as (4).
  6. Run the model with SubDyn and SoilDyn. Switch from SoilDyn CalcOption 1 to 2 and define a nonlinear stiffness.

I hope that helps.

Roger

1 Like

Dear Roger,

Thank you for the inputs.
I’ve been testing from the simplest scenario and adding complexity to the system.
However, I cannot move from the point number (2). When I include the SubDyn with a clamp at the base and releasing all the DoF’s in the platform (ElastoDyn) I get the error as shown below. The only way to move forward is to restrict some of the DoF’s in the platform, but the outcome of the simulation is not realistic, as shown before.

The SubDyn input file includes the following:

-------------------------- SIMULATION CONTROL  ---------------------------------
False            Echo        - Echo input data to "<rootname>.SD.ech" (flag)
"DEFAULT"        SDdeltaT    - Local Integration Step. If "default", the glue-code integration step will be used.
             3   IntMethod   - Integration Method [1/2/3/4 = RK4/AB4/ABM4/AM2].
True             SttcSolve   - Solve dynamics about static equilibrium point
True             GuyanLoadCorrection - Include extra moment from lever arm at interface in interface reactions. R. Bergua: This accounts for the gravity being applied over the deflected structure.
-------------------- FEA and CRAIG-BAMPTON PARAMETERS---------------------------
             3   FEMMod      - FEM switch: element model in the FEM. [1= Euler-Bernoulli(E-B);  2=Tapered E-B (unavailable);  3= 2-node Timoshenko;  4= 2-node tapered Timoshenko (unavailable)]
             1   NDiv        - Number of sub-elements per member
True             CBMod       - [T/F] If True perform C-B reduction, else full FEM dofs will be retained. If True, select Nmodes to retain in C-B reduced system.
             0   Nmodes      - Number of internal modes to retain (ignored if CBMod=False). If Nmodes=0 --> Guyan Reduction.
			 1   JDampings   - Damping Ratios for each retained mode (% of critical) If Nmodes>0, list Nmodes structural damping ratios for each retained mode (% of critical), or a single damping ratio to be applied to all retained modes. (last entered value will be used for all remaining modes).
             0   GuyanDampMod - Guyan damping [0=none, 1=Rayleigh Damping, 2= user specified 6x6 matrix]
0, 0             RayleighDamp - Mass and stiffness proportional damping  coefficients (Rayleigh Damping) [only if GuyanDampMod=1]
             6   GuyanDampSize  - Guyan damping matrix (6x6) [only if GuyanDampMod=2]. Roger comment: Damping matrix applied at the Transition Piece. 
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
---- STRUCTURE JOINTS: joints connect structure members (~Hydrodyn Input File)---
            2   NJoints     - Number of joints (-). Roger comments: JointType: (1): cantilever beam, (2) universal joint, (3) revolute joint, (4) spherical joint. 
JointID          JointXss               JointYss               JointZss          JointType     JointDirX   JointDirY     JointDirZ    JointStiff    JointDamp   ![Coordinates of Member joints in SS-Coordinate System]
  (-)               (m)                    (m)                    (m)               (-)           (-)         (-)           (-)        (Nm/rad)     (N/rad.s)
   1                 0                      0                      0               1			   0		   0			 0			   0			0
   2                 0                      0                     -0.2             1		   	   0	       0		     0		       0  		    0    ! Revolute joint        
------------------- BASE REACTION JOINTS: 1/0 for Locked/Free DOF @ each Reaction Node ---------------------
             1   NReact      - Number of Joints with reaction forces; be sure to remove all rigid motion DOFs of the structure  (else det([K])=[0])
RJointID   RctTDXss    RctTDYss    RctTDZss    RctRDXss    RctRDYss    RctRDZss   SSIfile ![Global Coordinate System] 
  (-)       (flag)      (flag)      (flag)      (flag)      (flag)      (flag)    (string)
   2          1           1           1            1          1            1  
------- INTERFACE JOINTS: 1/0 for Locked (to the TP)/Free DOF @each Interface Joint (only Locked-to-TP implemented thus far (=rigid TP)) ---------
             1   NInterf     - Number of interface joints locked to the Transition Piece (TP):  be sure to remove all rigid motion dofs
IJointID   ItfTDXss    ItfTDYss    ItfTDZss    ItfRDXss    ItfRDYss    ItfRDZss   
  (-)       (flag)      (flag)      (flag)      (flag)      (flag)      (flag)
   1           1           1           1           1           1           1
----------------------------------- MEMBERS --------------------------------------
              1   NMembers    - Number of frame members. 
MemberID   MJointID1   MJointID2   MPropSetID1   MPropSetID2     MType        COSMID        Roger comments: MType: (1): beam, (2) pre-tensioned cable, (3) rigid connection.
  (-)         (-)         (-)          (-)           (-)          (-)              (-)
   1           1           2            1             1            1                			          
------------------ MEMBER X-SECTION PROPERTY data 1/2 [isotropic material for now: use this table for circular-tubular elements] ------------------------
              1   NPropSets   - Number of structurally unique x-sections
PropSetID     YoungE          ShearG          MatDens          XsecD           XsecT
  (-)         (N/m2)          (N/m2)          (kg/m3)           (m)             (m)
   1  		  4.2E10		  1.75E10		   1375			    5.1		         0

Considering only the ElastoDyn with the DoF’s fixed at the platform I get consistent results without any error, as shown by the plot below:
2

In my view, for some reason, SubDyn is not connecting properly with ElastoDyn at the interface joint.

Thank you.

Best regards,
Mário

Dear @Mario.Sobrinho,

It’s good to know that the issue is coming from the beam introduced in SubDyn.

I’m wondering if you are using a model available in a public GitHub repository. If so, I could try adding the beam in SubDyn and figuring out what is going on.

Clearly, your model goes numerically unstable when including the small beam in SubDyn. Some things that you can try to avoid this instability are:

  1. Reduce the time step in OpenFAST.
  2. Include 1 correction iteration in OpenFAST input file.
  3. Include a damping for the beam that you are introducing. For example: RayleighDamp = 0, 0.0106. This would be a Rayleigh stiffness proportional only damping providing 1% critical damping at the frequency of 0.3 Hz and growing linearly for higher frequencies.
  4. In SubDyn note that XsecT is the beam thickness. I understand that you are trying to define a solid beam?
  5. The beam introduced in SubDyn seems to be very short (0.2 m) and wide (5.1 m in diameter). I expect this to result in a very high frequency for this SubDyn part (this frequency can be checked in the SubDyn summary file, *.yaml extension: look for “Frequencies of Guyan modes”). This could challenge the ElastoDyn-SubDyn coupling requiring quite small time steps. Actually, a message should be printed in the console when running OpenFAST about this frequency and a recommended time step.
  6. I do not expect this short and stiff beam to change the dynamics in your system much. However, note that if the dynamics change (e.g., introducing some compliance at the bottom like a stiffness matrix at the seabed), you should compute again the mode shapes for the tower according to the new boundary conditions. ElastoDyn relies in a assumed modes approach where the user must provide the tower mode shapes as input.

I hope that helps!

Roger

Dear Roger,

Thank you for the inputs.
I believe I was able to understand what I was doing wrong, by reading your point number 4. I was not paying attention to the MEMBER SECTION PROPERTY. For now, by changing the dimensions of XsecD and XsecT the model was able to run, even considering the REDWIN DLL model 1, with results that seems to be in line with the expected.
Once again, thank you for your precious inputs. I hope not need to bother you again with this topic :slight_smile:

Thank you and best regards,
Mário

2 Likes