Can Elastodyn directly consider hydrodynamic loads

Dear all,

I am confused whether Elastodyn can directly consider hydrodynamic loads?

For offshore turbine fixed by monopile, I know we can use Subdyn to model monopile and tower modelled in Elastodyn can be mounted on monoplie. And Subdyn can consider hydrodynamic loads through mapping.

But I guess I can use Elastodyn to model the combination of tower and monopile if I know the combination’s mode shape. But I am not sure whether the hydrodynamic loads in Hydrodyn can be correctly mapped to Elastodyn.

The reason I want to only use Elastodyn is because:

  1. Elastodyn is more efficient than Subdyn.
  2. The damping definitions are different between Elasotdyn and Subdyn. Elastodyn adopts critical damping ratio. While Subdyn uses a combination of CB damping and Rayleigh damping. I want to avoid using different damping.


Dear @Ran.Tu,

Neither FAST v8 nor OpenFAST support a mesh-mapping between the tower of ElastoDyn and HydoDyn. So, without changing the source code, if you want the hydrodynamic loads to be applied to a flexible substructure, you must enable and represent the substructure in SubDyn.

Best regards,

Thank you Jason,

Now I am simulating offshore fixed wind turbine supporting by monopile. When I use Subdyn to simulate monopile, the time increment in .fst file should decrease to 0.000625 to avoid divergence. For all modules (including Aerodyn15, servodyn, elastodyn), I used the default time increment (identical to 0.000625). It took a lot of time to complete the simulation as such small time increment was adopted.

  1. So I am wondering if there are some skills which could speed up the calculation?

  2. What is the linear algebra solver adopted by OpenFAST? I know there are some open source linear algebra packages, i.e. lapack, Eigen, PETSc and Trilinos, which are powerful and can be introduced into OpenFAST. It seems that Map++ has adopted lapack. But I am not sure whether other modules also adopt some efficient linear algebra solvers.
    update: After roughly looking at Subdyn source code, I found Subdyn also adopts lapack


Dear @Ran.Tu,

The time step you are using seems quite small. How did you determine that this time step is needed? How many Craig-Bampton modes have you enabled in SubDyn and what is the highest Craig-Bampton mode frequency? Some guidance for setting up an offshore model in OpenFAST with SubDyn enabled are given here:

Regarding linear algebra solvers, this likely depends on which module you are asking about, but I agree that LAPACK is likely the most frequently used.

Best regards,

Dear Jason:

I modelled a 30m long monopile in Subdyn. The outer diameter of the monopile is 6m, and the thickness is 0.06m. The monopile is splited into 3 members. And each member is further divided into 3 elements. Timoshenko beam is adopted. Below is my damping:

             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)]
             3   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.
             4   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).
             1   GuyanDampMod - Guyan damping {0=none, 1=Rayleigh Damping, 2=user specified 6x6 matrix}
  0.000, 0.010   RayleighDamp - Mass and stiffness proportional damping  coefficients (Rayleigh Damping) [only if GuyanDampMod=1]
             6   GuyanDampSize - Guyan damping matrix (6x6) [only if GuyanDampMod=2]
   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   1.4625e+06   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

After running simulation, I got below information:

Running SubDyn.
   Fixed bottom case detected
   Performing Craig-Bampton reduction 60 DOFs -> 4 modes + 6 DOFs
   Using static improvement method for gravity and ext. loads
   SubDyn recommended dt: 0.00090 - Current dt: 0.00063 - Max frequency:   55.764
   Calculating Full System Modes for output files
   Exporting Summary file

As you can see, the recommended dt = 0.00090. So I choose 0.000625 to avoid divergency. Below is some subdyn infomation. Hope it can reflect something…

#Eigenfrequencies [Hz] for full system, with reaction constraints (+ Soil K/M + SoilDyn K0) 
Full_frequencies: # 1 x 54
  - [   6.319104E+00,   6.319104E+00,   2.676944E+01,   2.946952E+01,   2.946952E+01,   4.315620E+01,   6.096864E+01,   6.096864E+01,   8.112554E+01,   8.678901E+01,   8.678901E+01,   1.060086E+02,   1.060086E+02,   1.191718E+02,   1.191718E+02,   1.279464E+02,   1.279464E+02,   1.307861E+02,   1.334173E+02,   1.334173E+02,   1.364113E+02,   1.364113E+02,   1.379490E+02,   1.988829E+02,   2.223937E+02,   2.653209E+02,   3.206279E+02,   3.375798E+02,   3.756434E+02,   3.756434E+02,   4.128127E+02,   4.277356E+02,   4.380972E+02,   4.380972E+02,   4.813307E+02,   5.227480E+02,   5.227480E+02,   5.246536E+02,   5.442274E+02,   6.008514E+02,   6.008514E+02,   6.628179E+02,   6.628179E+02,   6.655137E+02,   7.081735E+02,   7.081735E+02,   7.391388E+02,   7.391388E+02,   7.578514E+02,   7.578514E+02,   7.759746E+02,   7.764889E+02,   7.764889E+02,   8.458174E+02,]
#Frequencies of Guyan modes [Hz]
GY_frequencies: # 1 x 6
  - [   6.391218E+00,   6.391218E+00,   2.948010E+01,   4.752618E+01,   5.409784E+01,   5.409784E+01,]
#Frequencies of Craig-Bampton modes [Hz]
CB_frequencies: # 1 x 4
  - [   2.778162E+01,   2.778162E+01,   5.374285E+01,   5.576394E+01,]


In Subdyn theory, I found:

Until full-system linearization is made available, experience has shown that it is sufficient to enable all C-B modes up to 10 Hz (the natural frequencies of the C-B modes are written to the SubDyn summary file)

I do not fully understand. But for my model, the lowest frequency of C-B frequency is 27.78 Hz which is larger than 10 Hz. Does it mean there is no need to retain 4 modes? But how many modes should I retain?


Dear @Ran.Tu,

Our recommended rule of thumb for choosing the time step is DT = 1/(10 * fmax), where fmax is the highest natural frequency of the model in Hz. In your case, this would be DT = 1/10*55.76 Hz ) = 0.0017 s. (I’m not fully sure why SubDyn recommends using half this value). Regardless, I would think 0.000625 s would be plenty small, assuming you don’t have much higher frequencies in ElastoDyn. Have you applied 1 correction step (NumCrctn = 1 in the OpenFAST primary input file), which is often required when enabling SubDyn per the modeling guidance I linked to above?

Regardless, I agree with the SubDyn documentation recommendation of keeping only the lowest-frequency Craig-Bampton, less than 10 Hz. In your case, it means you should be able get by without any Craig-Bampton modes, Nmodes = 0 in SubDyn, which means that you’ll likely be able to increase the time step quite a bit.

Best regards,

Thank you Jason.

  1. For NumCrctn, I have tried 0 and 1. All the results except reaction force in z axis does not change. Does NumCrctn refer to SIM? After reading Subdyn theory, I found SIM ( Static-Improvement Method) will influence the calculation of gravity and bouyancy.

I checked the Subdyn modeling consideration. And I found for the RK4 and ABM4 methods, DT = 1/(10 * fmax). But for the AB4 method, the recommended time step is half DT. Weird thing is I have chose ABM4. As you can see below:

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
False            GuyanLoadCorrection - Include extra moment from lever arm at interface and rotate FEM for floating.

After reading more of Subdyn modeling consideration, I get some questions:
3. It said:
however, for tubular towers, the structural model in ElastoDyn tends to be more accurate because ElastoDyn considers geometric nonlinearities not treated in SubDyn.

I am confused by “geometric nonlinearities” of Elastodyn. Here do you mean P-delta effect? In my understanding, P-delta effect will lead to extra moment because of the offset between undeflected position and deflected position. But according to subdyn theory, P-delta effect can be accounted by enabling GuyanLoadCorrection.

  1. It said:
    If FAST is run with SubDyn but not HydroDyn, the water depth will be automatically set to 0 m. This will influence the calculation of the reaction loads. Reactions are always provided at the assumed mudline, therefore, they would not be correctly located for an offshore turbine as a result.

My understanding is still water level is always located at z = 0. If Hydrodyn is turned off, the water depth will be automatically set to 0 m. That is to say, mudline is located at z = 0 ( the same with still water level). And OpenFAST always assumes reaction force located at mudline (z = 0). But if I set reaction point at z=-20 in Subdyn file, this discrepancy (I manually set z of reaction point = -20, but OpenFAST assume z of reaction point = 0) will lead to incorrect calculation.

In Subdyn summary file, I get Mbbt and Kbbt. Does Mbbt and Kbbt refer to MBB and KBB respectively shown in above figure?

Der @Ran.Tu,

Here are my responses:

  1. No, I was referring to the number of correction steps (NumCrctn) specified in the OpenFAST primary input file, which is independent of the static improvement method (SIM) within SubDyn. The OpenFAST glue code implements a loose-coupling approach based on a predictor-corrector formulation, with NumCrctn determining the number of correction steps. For most models with SubDyn enabled, we typically recommend setting NumCrctn = 1.
  2. Looking at the SubDyn source code, I see that the “SubDyn recommended dt” is calculated as 1/(20 * freq_max), regardless of the integrator specified. I would have expected this to be 1/(10 * freq_max). Regardless, it is in the ballpark.
  3. The P-delta effect is one term. The other term is the “axial shortening” of the tower–that is, the downwards deflection of the tower associated with transverse deflection because the tower does not stretch when deflected transversely. This effect is important for gravitational destiffening associated with large tower-top mass. (For blades, this effect is important for centrifugal stiffening). This term is accounted for in ElastoDyn, but not SubDyn. Despite ElastoDyn using mode shapes as shape functions, ElastoDyn does consider several geometric nonlinear terms that improve its accuracy for small and moderate deflections relative to what a purely linear model would calculate.
  4. This statement from the SubDyn manual is just a warning regarding how to interpret the reaction load write outputs from SubDyn (ReactMxss, etc.) when HydroDyn is disabled for an offshore wind turbine. Please note that as of OpenFAST v3.1 and newer, however, the WtrDpth is now specified in the OpenFAST primary input file, and so, can be set nonzero even when HydroDyn is disabled. So, I’m not sure this statement from the SubDyn manual is really valid anymore for newer versions of OpenFAST.
  5. Yes, except that the “t” in the names (Mbbt and Kbbt) refer to the versions of M_BB and K_BB with a tilde, i.e., after the boundary constraints are applied (see: SubDyn Theory — OpenFAST v3.2.0 documentation).

Best regards,

Thank you Jason for the explanation,

Now I am studying the subdyn source code. And I arise two questions during reading the code.

  1. In FEM.f90, subroutine ElemK_Beam and ElemM_Beam calculate beam element mass matrix and beam element stiffness matrix respectively. They both depends on directional cosine matrix. During simulation, elements will diplace/deform, and can change directional cosine matrix, leading to a change of both beam element mass matrix and beam element stiffness matrix. So I am curious that will beam element mass matrix and beam element stiffness matrix be updated during simulation? Even for the lumped forces and moments due to gravity, will it be updated?

  2. In SD_FEM.f90, I found Subdyn can consider both SSI matrix and soil matrix. I know SSI matrix is inputted for reaction point. But I don’t find where soil matrix (it is called Init%Soil_K in code) is inputted in subdyn primary file.

Dear @Ran.Tu,

Regarding (1), what you describe is a geometric nonlinear effect (the deflection impacting the mass and stiffness matrices) not considered in SubDyn.

Regarding (2), Soil_K is an initialization input to SubDyn, i.e., an input defined by the calling program. When SubDyn is coupled to OpenFAST, Soil_K is set by OpenFAST (and is not currently set in OpenFAST v3.2). Once the SoilDyn pull request is merged in (SoilDyn -- soil dynamics module by andrew-platt · Pull Request #986 · OpenFAST/openfast · GitHub), Soil_K will be an initialization output of SoilDyn and transferred as an initialization input to SubDyn through the OpenFAST glue code.

Best regards,