Restoring Force upon Steady state in WEC-Sim

[Elaine.Liu]

Hi all,
I am running a simulation of an offshore WEC with its mooring line, analyzed using MoorDyn. I am trying to implement a delayed wave ramp function that applies waves onto the structure after a certain time has elapsed, i.e. 120 seconds. I expect that there is some free decay in the beginning, and that it would stay at its natural position (given that there is no wave excitation), yet there are small oscillations that arise. It seems these contributions are from the restoring/mooring force. This does not damp out over time as well. Is there a way to reach full steady state of the object, where no motions are occurring, and then apply the wave?

Thanks for the help!

Best,
Elaine

[Mohamed.Shabara]

It looks like there might be an issue with the hydrodynamic coefficients generated by your BEM solver. In the free decay test, the motion does not decay to zero, which suggests a potential problem with the hydrodynamic stiffness or radiation damping coefficients.

I recommend checking the following:

1. Mesh Quality: Ensure your mesh has low skewness and is of good quality. A highly skewed mesh can lead to inaccurate results.
2. Hydrodynamic Coefficient Convergence: Use BEMIO to check whether your added mass and radiation damping plots converge to constant values. If they do not, try increasing the maximum frequency input in your BEM solver until convergence is achieved.
3. Free Decay Test: Run a separate free decay test using: wave(1) = waves('noWaveCIC'); before testing your model in wavy environment.

[Elaine.Liu]

Hi Mohamed,
Thank you for the insightful comments. I ran a free decay test with ‘noWave’ and here are my results. It seems that it does in fact reach steady state at a neutral position as shown in the heave direction. This means the oscillating behavior after I delay the wave ramping is strange. Do you think this still has to do with the hydrodynamic coefficients? Thanks!

Best,
Elaine

[Elaine.Liu]

Hi,
Just following up on this to see if you had any insight! I’m happy to send over my simulink model if that would be easier to visualize.

Thanks,
Elaine

[Mohamed.Shabara]

I assume you didn’t add an initial displacement? If no initial displacement was applied, the platform shouldn’t exhibit a -0.12 m motion. This suggests a potential issue with the system setup.

Additionally, your stiffness values appear to be quite small or negative in all directions except the heave motion. It’s well known that BEM tends to produce low stiffness values for yaw motion, which could explain why your floater continues to yaw indefinitely, even in the absence of an applied torque.

If you are confident that your BEM code is functioning correctly, consider the following:

1. Adjust the Center of Gravity (C.G.) – Ensure that the C.G. is set correctly to match the actual weight distribution. This should prevent any unintended motion in the heave direction.
2. Modify the Constraints – Try replacing the 6DOF constraint in your Simulink model with a 3DOF constraint that allows only x, y, and z motions. This will help eliminate excessive yaw motion.
3. Verify the BEM Code – If the problem persists even with a 3DOF constraint, there may be an issue with the BEM results that needs to be investigated further.

Let me know if you need further clarification or assistance in debugging the issue!

[Elaine.Liu]

Hi Mohamed,
That all makes sense, I will go ahead and investigate your suggestions, thanks!

Elaine

[Chaozhi.Qiu]

Dear @Mohamed.Shabara ,
Hello!
I am also conducting the OC3 spar free decay test in WEC-Sim and comparing it with OpenFAST.

In OpenFAST, I disabled AeroDyn, MoorDyn and only left the spar floater. And set ExctnMod = 0, also set the strip-theory related variables to zeros to ensure only consider the potential-flow, and I also zero out the AddCLin, AddBLin, and AddBQuad.
In the WEC-Sim, I set the wave(1) = waves('noWaveCIC');. The heave-free decay test matched well in both codes. However, the pitch-free decay test shows some frequency difference, as indicated in the figure (set 5-degree initial pitch angle).

For your reference, I have checked the added mass matrix and hydrostatic matrix. The hydrostatic matrix is exactly the same, but the added mass matrices in OpenFAST and WEC-Sim are different, as:
OpenFAST:
7.759112E+06 0.000000E+00 0.000000E+00 0.000000E+00 -4.831406E+08 0.000000E+00
0.000000E+00 7.759089E+06 0.000000E+00 4.831395E+08 0.000000E+00 0.000000E+00
0.000000E+00 0.000000E+00 2.412549E+05 0.000000E+00 0.000000E+00 0.000000E+00
0.000000E+00 4.831413E+08 0.000000E+00 3.793609E+10 0.000000E+00 0.000000E+00
-4.831428E+08 0.000000E+00 0.000000E+00 0.000000E+00 3.793618E+10 0.000000E+00
0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 2.598276E-06
WEC-Sim:
-2.375979+07 0.000000E+00 0.000000E+00 0.000000E+00 -4.83140e+08 0.000000E+00
0.000000E+00 -2.375982e+07 0.000000E+00 483139490 0.000000E+00 0.000000E+00
0.000000E+00 0.000000E+00 -3.127765e+07 0.000000E+00 0.000000E+00 0.000000E+00
0 483141335 0 0 0 0
-483142770
0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00

I checked the source code of WEC-Sim, and the added mass matrix was ‘lumped’ by considering the numerical stability.

My questions are:

• Will the lumped added mass cause the simulation results to differ? As the rotational DOFs in WEC-Sim added mass matrix are zero.
• In the WEC-Sim document, the governing equation,
  
  Now only Fb(t) and Frad(t) are included; the WEC-Sim pitch change frequency is smaller than OpenFAST. Are there any terms being ignored?
• For the surge motion, why does it move towards the negative X-axis, and the neutral position is not zero? (I am confused about the physical phenomenon.)

Any comments are welcome, and thanks in advance.

Best regards,
Chaozhi Qiu

[Mohamed.Shabara]

I noticed that in your free decay test, the response does not appear to decay over time for either solver. The pitch motion also suggests that the system in OpenFAST is moving faster than in WEC-Sim.

Regarding the added mass, could you share the BEMIO plots? Are you using the same BEM output for both solvers? Also, I wanted to confirm whether you are aware of the difference between noWaveCIC and noWave.

I’m not fully familiar with how the free decay test is computed in OpenFAST, but the differences in the added mass results seem like something worth comparing against WEC-Sim. I plan to take a closer look at this later this week. In addition, I noticed that your WEC-Sim added mass results include negative terms—can you confirm if that is expected?

Yours,

Mohamed

[Chaozhi.Qiu]

Thank you very much for the quick reply.

Because I zeroed out the AddCLin, AddBLin, and AddBQuad in OpenFAST, and set the LinearDamping damping matrix to zeros in WEC-Sim. Thus, the free decay test didn’t decay over time (Because I found the decay period in WEC-Sim is different from OpenFAST, I zeroed them out to separate the effects to check the decay period).
And I use the same BEM outputs for both solvers. The Spar.1, Spar.3, and Spar.hst files came from OpenFAST’s baseline data (OpenFAST hydro data). Then I use BEMRosetta to generate .out file for the WEC-Sim. The figures below are the BEMIO plots.

I am aware of the difference between noWaveCIC and noWave. The noWaveCIC uses a convolution integral calculation to calculate radiation damping, and the noWave uses the constant damping matrix to calculate. Am I right?

Regarding the added mass in OpenFAST and WEC-Sim. I output them from the solvers directly without modification, so the WEC-Sim added mass does contain negative terms. Additionally, I attempted to set the added mass as OpenFAST in WEC-Sim manually, but the WEC-Sim crashed due to an error: “The simulation will be stopped. There may be a singularity in the solution. If*
not, try reducing the step size (either by reducing the fixed step size or by tightening the error tolerances)”

Looking forward to your reply.

Best,
Chaozhi Qiu

[Mohamed.Shabara]

I am aware of the difference between noWaveCIC and noWave. The noWaveCIC uses a convolution integral calculation to calculate radiation damping, and the noWave uses the constant damping matrix to calculate. Am I right?

Yes, this is correct, but there is another difference too. The `noWaveCIC uses the added mass at infinit frequency $m_\infty$ while noWave uses the added mass at the frequency the user specify in the wecSimInoutFile.m

Regarding the added mass in OpenFAST and WEC-Sim. I output them from the solvers directly without modification, so the WEC-Sim added mass does contain negative terms. Additionally, I attempted to set the added mass as OpenFAST in WEC-Sim manually, but the WEC-Sim crashed due to an error: “The simulation will be stopped. There may be a singularity in the solution. If*
not, try reducing the step size (either by reducing the fixed step size or by tightening the error tolerances)”

Added mass is a special hydrodynamic effect that makes time-domain modeling tricky. The challenge is that the added mass force is proportional to acceleration, but acceleration is what the solver is trying to compute. This creates an algebraic loop: the acceleration depends on the added mass force, but the added mass force itself depends on acceleration.

The most robust way to break this loop is to combine the added mass matrix with the rigid body’s mass and inertia on the left-hand side of the governing equations. This removes the acceleration dependence from the force summation, leaving a clean system of equations to solve. OpenFAST follows this approach, which is why it can directly use the unmodified added mass matrix.

However, Simscape Multibody (used by WEC-Sim) does not allow this formulation. In Simscape, a rigid body can only have one translational mass and one inertia tensor. It cannot have direction-dependent masses (e.g., a different effective mass in surge vs heave). Because of this limitation, WEC-Sim cannot directly combine the added mass with the body mass and inertia.

To work around this, WEC-Sim:

1. Lumps certain added mass terms into the body mass/inertia — rotational added mass terms are fully mapped into the body inertia tensor, while the surge, sway, and heave diagonal terms are lumped into the translational mass.

2. Reduces the remaining added mass force — the “excess” added mass terms are left on the right-hand side, but with reduced magnitude.

3. This creates modified mass and added mass matrices that stabilize the time-domain simulation in Simscape while still capturing the physics.

Because of this implementation difference, if you directly overwrite WEC-Sim’s added mass with OpenFAST values, you can create a singular mass matrix that Simscape cannot solve — which is why your simulation crashed with the “singularity” error.

For more details, see the WEC-Sim documentation on added mass handling:
https://wec-sim.github.io/WEC-Sim/dev/developer/advanced_features.html#added-mass

[Jason.Jonkman]

Just one correction regarding OpenFAST: OpenFAST does not combine the added mass with the mass matrix on the left-hand side. OpenFAST keeps the added mass in the hydrodynamic loads (on the right-hand side, dependent on acceleration) and solves the force-acceleration coupling implicitly by effectively solving constraint equations (the input-output solve) at each time step. With this approach, there is no limitation to the values in the added mass matrix. Moreover, there is not a single 6x6 added mass matrix, but rather, support for multiple bodies in the potential-flow solution and hydro-elastic added mass in the strip-theory solution.

Best regards,

[Chaozhi.Qiu]

Hi, @Jason.Jonkman

Thank you very much for the correction.

Could you please explain how the OpenFAST “solves the force-acceleration coupling implicitly by effectively solving constraint equations (the input-output solve) at each time step.”?

Besides, as I asked above, I compared the OpenFAST with WEC-Sim for the OC3 spar free decay. And I set ExctnMod = 0 and RdtnMod = 0, and zero out the AddCLin, AddBLin, and AddBQuad, also set the Cd in the strip-theory to zero to ‘disable’ the morison elements calculation. And the mass and inertia of tower and blades also set to nearly zeros. To achieve a static equilibrium at heave=0, I adjusted the platform mass (PtfmMass) from the default ~7.47e6 kg to 8.23e6 kg. Therefore, I assumed the motion of equations should be simplified as:

And I code a simple example to simulate the pitch free decay, and put the results from OpenFAST, code, and WEC-Sim together as below:

While the pitch frequency is different in the three codes, the surge amplitudes also have a large difference.
My questions are:

• With my specified settings, is my simplified equation of motion correct for what OpenFAST is solving? Or are there other residual terms that might remain active.
• For the free decay test, radiation and excitation will play important role?
• How does changing PtfmMass affect the coupled dynamics? While it corrected the static heave, could this mass change be the reason for the discrepancy in the pitch natural frequency and the resulting surge motion.

Thank you for your time and any insight you can provide.

Best regards

[Chaozhi.Qiu]

Thank you very much for the detailed explanation.
And I uploaded the figure under the answer of @Jason.Jonkman. In the example code, I use the full added mass matrix (A_inf), and the WEC-Sim still uses the lumped added mass matrix. The results from WEC-Sim and my code are almost identical. However, the pitch natural frequency still does not match the OpenFAST.

I can’t think of any other difference in settings now.

Any further suggestions on what to investigate would be greatly appreciated.

Best regards

[Jason.Jonkman]

Dear @Chaozhi.Qiu,

I can’t comment on WEC-Sim, but I’ll comment on OpenFAST:

• We have a paper under development that explains the original implicit loose coupling scheme of OpenFAST (including implicit force-acceleration coupling), as well as the new tight coupling scheme that will be introduced soon in OpenFAST v5. We intend to submit this paper when OpenFAST v5 is released in the next month or so.
• By setting RdtnMod = 0, both the added mass and the wave-radiation damping will be zeroed, so, A_inf in your equation will also be zeroed.
• Your equation misses the gravitational force from body weight, i.e. PtfmMass * Gravity, unless you’ve lumped this into F_buoy. You haven’t said what you are doing with other system masses like yaw bearing, nacelle, and hub, but these would be included here as well.
• By eliminating component masses like tower and blades, and adding mass to the platform, this would change the center of mass of the system. This will have an impact on the gravitational restoring in pitch and roll. You haven’t stated what terms is included in C_static, but normally hydrostatic restoring includes contributions from both body weight and buoyancy, where the former in OpenFAST is accounted for in the structural module (ElastoDyn, BeamDyn, and SubDyn, when enabled) and the later is accounted for in the hydrodynamic module (HydroDyn). The shifting of body masses around will impact the former in OpenFAST.

Best regards,

[Chaozhi.Qiu]

Dear @Jason.Jonkman ,

Thank you very much for the reply.

• I am looking forward to reading your paper.

• I notice in the WAMIT_CalcOutput of WAMIT.f90:

  m%F_PtfmAM     =   -matmul(p%HdroAdMsI, qdotdot) ! In h-frame
  do iBody = 1, p%NBody
     indxStart = (iBody-1)*6+1
     indxEnd   = indxStart+2
     call hiFrameTransform( h2i, u%PtfmRefY, m%F_PtfmAM(indxStart:indxEnd), tmpVec3, ErrStat2, ErrMsg2 )
     m%F_PtfmAM(indxStart:indxEnd) = tmpVec3
  
     indxStart = indxEnd+1
     indxEnd   = indxStart+2
     call hiFrameTransform( h2i, u%PtfmRefY, m%F_PtfmAM(indxStart:indxEnd), tmpVec3, ErrStat2, ErrMsg2 )
     m%F_PtfmAM(indxStart:indxEnd) = tmpVec3
  end do
  
     ! Compute outputs here:
  do iBody = 1, p%NBody
     indxStart = (iBody-1)*6
     
     DO I=1,3
        y%Mesh%Force(I,iBody)    = m%F_PtfmAM(indxStart+I)   + m%F_Rdtn(indxStart+I)   + m%F_Waves1(indxStart+I)   + m%F_HS(indxStart+I) 
     END DO
     DO I=1,3
        y%Mesh%Moment(I,iBody)   = m%F_PtfmAM(indxStart+I+3) + m%F_Rdtn(indxStart+I+3) + m%F_Waves1(indxStart+I+3) + m%F_HS(indxStart+I+3)
     END DO
  end do
  

This code section, which is below the code sets the total load contribution from radiation damping to zero; it seems the added mass is included in the force computation. This is a little bit confusing.

• Sorry for the confusion about the equation. I have lumped gravitational force into F_buoy. But I didn’t add other system masses.
• Thanks for the clarification. And regarding the C_static, I use the data from OpenFAST, which only contains C33, C44, and C55. If the center of mass changes, do the C_static need to be transformed accordingly?

Thank you for your time and any insight you can provide.

Best regards,

[Jason.Jonkman]

Dear @Chaozhi.Qiu,

Yes, in OpenFAST, the added mass is included in the hydrodynamic loads (not moved/added to the body mass matrix), as I mentioned above.

Yes, C_static should be influenced by a change to the center of mass because this should include the contributions of both buoyancy and body weight.

Best regards,

[Chaozhi.Qiu]

Dear @Jason.Jonkman ,
Thank you very much for your kind reply.

I have some confusion about the WAMIT input file used in OpenFAST, which considers (0,0,0) as the body reference point (same as when the input file sets PtfmRefxt, PtfmRefyt, and PtfmRefyt to zeros). Right?
If so, my questions are:

• In my case, I set PtfmPitch to 5 degrees. The platform rotates about which point? CoG or WAMIT reference point (0,0,0)?
• If the rotation point is around CoG, need to transfer the hydrodynamics coefficients to reference at CoG? or the HydroDyn only calculate the hydro loads at (0,0,0), and the loads aggregation happens in ElastoDyn?

I hope I stated my questions clearly

Thank you for your time and any insight you can provide.

Best regards,

[Jason.Jonkman]

Dear @Chaozhi.Qiu,

The reference point of the WAMIT solution should be aligned with HydroDyn inputs PtfmRefxt, PtfmRefyt, and PtfmRefzt, which it sounds like should be (0,0,0) in your case.

PtfmPitch is an ElastoDyn input and is specified about ElastoDyn’s platform reference point, which is located by ElastoDyn input PtfmRefzt, which need not match HydroDyn’s potential flow reference point.

If the ElastoDyn reference point and HydroDyn potential flow reference point are not aligned, the coupling between ElastoDyn and HydroDyn will properly represent the kinematics, i.e., pitch motion about the ElastoDyn reference point will lead to both pitching and surging of the HydroDyn potental-flow reference point, with the moment arm equal to the offset between the two points. This is handled by the spatial mesh-to-mesh mapping between modules in the OpenFAST glue code.

Best regards,

[Chaozhi.Qiu]

Dear @Jason.Jonkman,

Thank you very much for the reply; it solves my confusion.

Best regards,