Platform submitted to regular waves : heave motion

Hello,

I currently work on a modeling of an offshore platform and I want to get the RAO of its heave motion when it’s only submitted to waves excitation.
I modelized my system with the potential flow theory, I put WaveMod on 1 because I want regular waves and the wave Hs is 2. I don’t use any second order, there is no current, no ballast and no marine growth. I used Wamit files (.1, .3 and .hst) in entry of Hydrodyn.
Here is my file Hydrodyn :

------- HydroDyn v2.03.* Input File --------------------------------------------
NREL 5.0 MW offshore baseline floating platform HydroDyn input properties for the OC4 Semi-submersible.
TRUE Echo - Echo the input file data (flag)
---------------------- ENVIRONMENTAL CONDITIONS --------------------------------
1025 WtrDens - Water density (kg/m^3)
200 WtrDpth - Water depth (meters)
0 MSL2SWL - Offset between still-water level and mean sea level (meters) [positive upward; unused when WaveMod = 6; must be zero if PotMod=1 or 2]
---------------------- WAVES ---------------------------------------------------
1 WaveMod - Incident wave kinematics model {0: none=still water, 1: regular (periodic), 1P#: regular with user-specified phase, 2: JONSWAP/Pierson-Moskowitz spectrum (irregular), 3: White noise spectrum (irregular), 4: user-defined spectrum from routine UserWaveSpctrm (irregular), 5: Externally generated wave-elevation time series, 6: Externally generated full wave-kinematics time series [option 6 is invalid for PotMod/=0]} (switch)
0 WaveStMod - Model for stretching incident wave kinematics to instantaneous free surface {0: none=no stretching, 1: vertical stretching, 2: extrapolation stretching, 3: Wheeler stretching} (switch) [unused when WaveMod=0 or when PotMod/=0]
3000 WaveTMax - Analysis time for incident wave calculations (sec) [unused when WaveMod=0; determines WaveDOmega=2Pi/WaveTMax in the IFFT]
0.05 WaveDT - Time step for incident wave calculations (sec) [unused when WaveMod=0; 0.1<=WaveDT<=1.0 recommended; determines WaveOmegaMax=Pi/WaveDT in the IFFT]
2 WaveHs - Significant wave height of incident waves (meters) [used only when WaveMod=1, 2, or 3]
251.3274 WaveTp - Peak-spectral period of incident waves (sec) [used only when WaveMod=1 or 2]
“DEFAULT” WavePkShp - Peak-shape parameter of incident wave spectrum (-) or DEFAULT (string) [used only when WaveMod=2; use 1.0 for Pierson-Moskowitz]
0.314159 WvLowCOff - Low cut-off frequency or lower frequency limit of the wave spectrum beyond which the wave spectrum is zeroed (rad/s) [unused when WaveMod=0, 1, or 6]
1.570796 WvHiCOff - High cut-off frequency or upper frequency limit of the wave spectrum beyond which the wave spectrum is zeroed (rad/s) [unused when WaveMod=0, 1, or 6]
0 WaveDir - Incident wave propagation heading direction (degrees) [unused when WaveMod=0 or 6]
0 WaveDirMod - Directional spreading function {0: none, 1: COS2S} (-) [only used when WaveMod=2,3, or 4]
1 WaveDirSpread - Wave direction spreading coefficient ( > 0 ) (-) [only used when WaveMod=2,3, or 4 and WaveDirMod=1]
1 WaveNDir - Number of wave directions (-) [only used when WaveMod=2,3, or 4 and WaveDirMod=1; odd number only]
0 WaveDirRange - Range of wave directions (full range: WaveDir +/- 1/2*WaveDirRange) (degrees) [only used when WaveMod=2,3,or 4 and WaveDirMod=1]
123456789 WaveSeed(1) - First random seed of incident waves [-2147483648 to 2147483647] (-) [unused when WaveMod=0, 5, or 6]
1011121314 WaveSeed(2) - Second random seed of incident waves [-2147483648 to 2147483647] (-) [unused when WaveMod=0, 5, or 6]
FALSE WaveNDAmp - Flag for normally distributed amplitudes (flag) [only used when WaveMod=2, 3, or 4]
“” WvKinFile - Root name of externally generated wave data file(s) (quoted string) [used only when WaveMod=5 or 6]
1 NWaveElev - Number of points where the incident wave elevations can be computed (-) [maximum of 9 output locations]
0 WaveElevxi - List of xi-coordinates for points where the incident wave elevations can be output (meters) [NWaveElev points, separated by commas or white space; usused if NWaveElev = 0]
0 WaveElevyi - List of yi-coordinates for points where the incident wave elevations can be output (meters) [NWaveElev points, separated by commas or white space; usused if NWaveElev = 0]
---------------------- 2ND-ORDER WAVES ----------------------------------------- [unused with WaveMod=0 or 6]
FALSE WvDiffQTF - Full difference-frequency 2nd-order wave kinematics (flag)
FALSE WvSumQTF - Full summation-frequency 2nd-order wave kinematics (flag)
0 WvLowCOffD - Low frequency cutoff used in the difference-frequencies (rad/s) [Only used with a difference-frequency method]
1.256637 WvHiCOffD - High frequency cutoff used in the difference-frequencies (rad/s) [Only used with a difference-frequency method]
0.618319 WvLowCOffS - Low frequency cutoff used in the summation-frequencies (rad/s) [Only used with a summation-frequency method]
3.141593 WvHiCOffS - High frequency cutoff used in the summation-frequencies (rad/s) [Only used with a summation-frequency method]
---------------------- CURRENT ------------------------------------------------- [unused with WaveMod=6]
0 CurrMod - Current profile model {0: none=no current, 1: standard, 2: user-defined from routine UserCurrent} (switch)
0 CurrSSV0 - Sub-surface current velocity at still water level (m/s) [used only when CurrMod=1]
“DEFAULT” CurrSSDir - Sub-surface current heading direction (degrees) or DEFAULT (string) [used only when CurrMod=1]
0 CurrNSRef - Near-surface current reference depth (meters) [used only when CurrMod=1]
0 CurrNSV0 - Near-surface current velocity at still water level (m/s) [used only when CurrMod=1]
0 CurrNSDir - Near-surface current heading direction (degrees) [used only when CurrMod=1]
0 CurrDIV - Depth-independent current velocity (m/s) [used only when CurrMod=1]
0 CurrDIDir - Depth-independent current heading direction (degrees) [used only when CurrMod=1]
---------------------- FLOATING PLATFORM --------------------------------------- [unused with WaveMod=6]
1 PotMod - Potential-flow model {0: none=no potential flow, 1: frequency-to-time-domain transforms based on WAMIT output, 2: fluid-impulse theory (FIT)} (switch)
“HydroDatabis\HydroData” PotFile - Root name of potential-flow model data; WAMIT output files containing the linear, nondimensionalized, hydrostatic restoring matrix (.hst), frequency-dependent hydrodynamic added mass matrix and damping matrix (.1), and frequency- and direction-dependent wave excitation force vector per unit wave amplitude (.3) (quoted string) [MAKE SURE THE FREQUENCIES INHERENT IN THESE WAMIT FILES SPAN THE PHYSICALLY-SIGNIFICANT RANGE OF FREQUENCIES FOR THE GIVEN PLATFORM; THEY MUST CONTAIN THE ZERO- AND INFINITE-FREQUENCY LIMITS!]
1 WAMITULEN - Characteristic body length scale used to redimensionalize WAMIT output (meters) [only used when PotMod=1]
11477 PtfmVol0 - Displaced volume of water when the platform is in its undisplaced position (m^3) [only used when PotMod=1; USE THE SAME VALUE COMPUTED BY WAMIT AS OUTPUT IN THE .OUT FILE!]
0 PtfmCOBxt - The xt offset of the center of buoyancy (COB) from the platform reference point (meters) [only used when PotMod=1]
0 PtfmCOByt - The yt offset of the center of buoyancy (COB) from the platform reference point (meters) [only used when PotMod=1]
1 ExctnMod - Wave Excitation model {0: None, 1: DFT, 2: state-space} (switch) [only used when PotMod=1; STATE-SPACE REQUIRES *.ssexctn INPUT FILE]
1 RdtnMod - Radiation memory-effect model {0: no memory-effect calculation, 1: convolution, 2: state-space} (switch) [only used when PotMod=1; STATE-SPACE REQUIRES *.ss INPUT FILE]
60 RdtnTMax - Analysis time for wave radiation kernel calculations (sec) [only used when PotMod=1 and RdtnMod>0; determines RdtnDOmega=Pi/RdtnTMax in the cosine transform; MAKE SURE THIS IS LONG ENOUGH FOR THE RADIATION IMPULSE RESPONSE FUNCTIONS TO DECAY TO NEAR-ZERO FOR THE GIVEN PLATFORM!]
“DEFAULT” RdtnDT - Time step for wave radiation kernel calculations (sec) [only used when PotMod=1 and RdtnMod=1; DT<=RdtnDT<=0.1 recommended; determines RdtnOmegaMax=Pi/RdtnDT in the cosine transform]
---------------------- 2ND-ORDER FLOATING PLATFORM FORCES ---------------------- [unused with WaveMod=0 or 6, or PotMod=0 or 2]
0 MnDrift - Mean-drift 2nd-order forces computed {0: None; [7, 8, 9, 10, 11, or 12]: WAMIT file to use} [Only one of MnDrift, NewmanApp, or DiffQTF can be non-zero]
0 NewmanApp - Mean- and slow-drift 2nd-order forces computed with Newman’s approximation {0: None; [7, 8, 9, 10, 11, or 12]: WAMIT file to use} [Only one of MnDrift, NewmanApp, or DiffQTF can be non-zero. Used only when WaveDirMod=0]
0 DiffQTF - Full difference-frequency 2nd-order forces computed with full QTF {0: None; [10, 11, or 12]: WAMIT file to use} [Only one of MnDrift, NewmanApp, or DiffQTF can be non-zero]
0 SumQTF - Full summation -frequency 2nd-order forces computed with full QTF {0: None; [10, 11, or 12]: WAMIT file to use}
---------------------- FLOATING PLATFORM FORCE FLAGS -------------------------- [unused with WaveMod=6]
True PtfmSgF - Platform horizontal surge translation force (flag) or DEFAULT
True PtfmSwF - Platform horizontal sway translation force (flag) or DEFAULT
True PtfmHvF - Platform vertical heave translation force (flag) or DEFAULT
True PtfmRF - Platform roll tilt rotation force (flag) or DEFAULT
True PtfmPF - Platform pitch tilt rotation force (flag) or DEFAULT
True PtfmYF - Platform yaw rotation force (flag) or DEFAULT
---------------------- PLATFORM ADDITIONAL STIFFNESS AND DAMPING --------------
0 0 0 0 0 0 AddF0 - Additional preload (N, N-m)
0 0 0 0 0 0 AddCLin - Additional linear stiffness (N/m, N/rad, N-m/m, N-m/rad)
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 AddBLin - Additional linear damping(N/(m/s), N/(rad/s), N-m/(m/s), N-m/(rad/s))
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 AddBQuad - Additional quadratic drag(N/(m/s)^2, N/(rad/s)^2, N-m(m/s)^2, N-m/(rad/s)^2)
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
---------------------- AXIAL COEFFICIENTS --------------------------------------
0 NAxCoef - Number of axial coefficients (-)
AxCoefID AxCd AxCa AxCp
(-) (-) (-) (-)
---------------------- MEMBER JOINTS -------------------------------------------
0 NJoints - Number of joints (-) [must be exactly 0 or at least 2]
JointID Jointxi Jointyi Jointzi JointAxID JointOvrlp [JointOvrlp= 0: do nothing at joint, 1: eliminate overlaps by calculating super member]
(-) (m) (m) (m) (-) (switch)
---------------------- MEMBER CROSS-SECTION PROPERTIES -------------------------
3 NPropSets - Number of member property sets (-)
PropSetID PropD PropThck
(-) (m) (m)
1 15.0000 0.0728 ! Base Column
2 10.0000 0.26 ! Column
3 5.2000 0.0282 ! Pontoons
---------------------- SIMPLE HYDRODYNAMIC COEFFICIENTS (model 1) --------------
SimplCd SimplCdMG SimplCa SimplCaMG SimplCp SimplCpMG SimplAxCa SimplAxCaMG SimplAxCp SimplAxCpMG
(-) (-) (-) (-) (-) (-) (-) (-) (-) (-)
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
---------------------- DEPTH-BASED HYDRODYNAMIC COEFFICIENTS (model 2) ---------
0 NCoefDpth - Number of depth-dependent coefficients (-)
Dpth DpthCd DpthCdMG DpthCa DpthCaMG DpthCp DpthCpMG DpthAxCa DpthAxCaMG DpthAxCp DpthAxCpMG
(m) (-) (-) (-) (-) (-) (-) (-) (-) (-) (-)
---------------------- MEMBER-BASED HYDRODYNAMIC COEFFICIENTS (model 3) --------
0 NCoefMembers - Number of member-based coefficients (-)
MemberID MemberCd1 MemberCd2 MemberCdMG1 MemberCdMG2 MemberCa1 MemberCa2 MemberCaMG1 MemberCaMG2 MemberCp1 MemberCp2 MemberCpMG1 MemberCpMG2 MemberAxCa1 MemberAxCa2 MemberAxCaMG1 MemberAxCaMG2 MemberAxCp1 MemberAxCp2 MemberAxCpMG1 MemberAxCpMG2
(-) (-) (-) (-) (-) (-) (-) (-) (-) (-) (-) (-) (-) (-) (-) (-) (-) (-) (-) (-) (-) ! Base of main column
-------------------- MEMBERS -------------------------------------------------
0 NMembers - Number of members (-)
MemberID MJointID1 MJointID2 MPropSetID1 MPropSetID2 MDivSize MCoefMod PropPot [MCoefMod=1: use simple coeff table, 2: use depth-based coeff table, 3: use member-based coeff table] [ PropPot/=0 if member is modeled with potential-flow theory]
(-) (-) (-) (-) (-) (m) (switch) (flag)
---------------------- FILLED MEMBERS ------------------------------------------
0 NFillGroups - Number of filled member groups (-) [If FillDens = DEFAULT, then FillDens = WtrDens; FillFSLoc is related to MSL2SWL]
FillNumM FillMList FillFSLoc FillDens
(-) (-) (m) (kg/m^3)
---------------------- MARINE GROWTH -------------------------------------------
0 NMGDepths - Number of marine-growth depths specified (-)
MGDpth MGThck MGDens
(m) (m) (kg/m^3)
---------------------- MEMBER OUTPUT LIST --------------------------------------
0 NMOutputs - Number of member outputs (-) [must be < 10]
MemberID NOutLoc NodeLocs [NOutLoc < 10; node locations are normalized distance from the start of the member, and must be >=0 and <= 1] [unused if NMOutputs=0]
(-) (-) (-)
---------------------- JOINT OUTPUT LIST ---------------------------------------
0 NJOutputs - Number of joint outputs [Must be < 10]
0 JOutLst - List of JointIDs which are to be output (-)[unused if NJOutputs=0]
---------------------- OUTPUT --------------------------------------------------
True HDSum - Output a summary file [flag]
False OutAll - Output all user-specified member and joint loads (only at each member end, not interior locations) [flag]
2 OutSwtch - Output requested channels to: [1=Hydrodyn.out, 2=GlueCode.out, 3=both files]
“ES11.4e2” OutFmt - Output format for numerical results (quoted string) [not checked for validity!]
“A11” OutSFmt - Output format for header strings (quoted string) [not checked for validity!]
---------------------- OUTPUT CHANNELS -----------------------------------------
“Wave1Elev”
“WRPSurge”
“WRPSway”
“WRPHeave”
“WRPRoll”
“WRPPitch”
“WRPYaw”
END of output channels and end of file. (the word “END” must appear in the first 3 columns of this line)

Concerning Elastodyn, I just enable the 6 degrees of freedom of platform motions, I put all the initial conditions to 0 and for the remaining parameters I just put the values according to my wind turbine and platform configuration.
When I launch a test, I just compute the Hydrodyn and Elastodyn files, and later I will add the moordyn module if I have solved the current problem.

At a very high period, I think I should have a heave motion with an amplitude of one, which is not the case here. For example when the wave period is about 251s, I get the “Heave_motion” file. I also created the .ss file of wamit in order to use a state space model for the radiation. I added this file to my wamit files and I put the RdtnMod on 2 instead of 1, but I still have the same problem with the “Heave_motion_ss” file.

It would be very helpful if you have an idea of the origin of this issue and how to correct it.

Best regards,

Julie Marchebois


Heave_motion.jpg

Dear Julie,

It sounds like you are running OpenFAST fine. You haven’t said anything about the floating platform you are simulating; why do you think the heave amplitude should be unity? Do you get the response amplitude that you expect at other frequencies?

Please note a couple things:

  • I see from your heave response that you have a mean heave displacement of about -0.5 m. This means that the buoyancy is not perfectly balanced with the weight of the floater, which is something I’d recommend that you fix.
  • It is likely easier to derive an RAO from OpenFAST by using white-noise wave excitation (WaveMod = 3), as discussed in our ISOPE 2013 paper: nrel.gov/docs/fy13osti/58098.pdf.

Best regards,

Dear Mr,

Thank you again for your quick answer.

-I was able to correct my problem of buoyancy and weight which were not well balanced: I have indeed made a mistake when I put values of some masses in Elastodyn. Now my heave motion is centered.

-Concerning my platform, I study a semi-submersible platform on which there is the tower of the wind turbine, but I don’t take into account the blades, the nacelle and the other elements in my model. I don’t know if I’m right but I put the masses/inertias of these elements to zero to eliminate their contribution.

-I thought that the response Amplitude Operators was the ratio of the vessel motion amplitude, here the heave amplitude, to incident wave amplitude which is equal to 1 in my case as the WaveHs is 2. This is why I thought that the heave motion was going to have an amplitude of 1 meter for low frequencies because I read that it’s the characterisitc behaviour of a heave motion : all the RAO of heave I saw tend to 1 for low frequencies, and the RAO of heave that I get thanks to the real experiments agrees with this observation.
I tried to compute my RAO anyway and for high frequencies I have a RAO which tends to 0 and I have a clear peak at a frequency which is, I guess, the natural frequency which seems to me quite coherent. However my graph is not as smooth as the graph on your document which is certainly due to the fact that I don’t have enough data for the moment. I’m also wondering if the total run time of my simulations is long enough (1500s in my FST file). And finally I still don’t really understand why a RAO is easier to obtain with a white noise spectrum of wave ? Can I still get a satisfying RAO by using regular waves on openfast ?

Best regards,
Julie Marchebois

Dear Julie,

OK, that all sounds reasonable. Does the RAO approach tend toward unity at lower and lower frequencies?

You can compute RAOs from OpenFAST by regular wave excitation, but this is not without challenges. One challenge associated with computing RAOs by regular waves is that OpenFAST is a nonlinear model, so, the response might not be sinusuidal and there may be response at frequencies different from the wave-excitation frequency. Another challenge is that it can be quite computationally expensive because you’ll need to run a long OpenFAST simulation at each wave frequency that you want to compute the RAO for. (You should be able to decide if the simulation is long enough by ensuring that it has reached a periodic solution.)

Computing RAOs by white noise is much more convenient because you can capture a large range of wave frequencies from a single simulation.

Best regards,

Dear Jason,

With regular waves, for lower and lower frequencies, it doesn’t tend to unity and I can’t find the origin of the problem.

Thank you for your advice, I now use the WaveMod 3 and I keep the same conditions than in my last simulations. I tried to follow step by step the method described in the document https://www.nrel.gov/docs/fy13osti/58098.pdf at the end of the page 3. I ran 6 simulations of 2500s with the WaveMod = 3 : for each simulation I put a random WaveSeed(1) and WaveSeed(2). I built two vectors, one containing the mean wave elevation of the 6 simulations (Wsum) and the other containing the mean heave displacement of the 6 simulations (Hsum) . I wrote a code to get the PSD of the Wsum and the CSD of Hsum and Wsum. Then I get the RAO by dividing my CSD by my PSD.
It seems that I made a mistake because the RAO obtained with white noise, as well as the RAO obtained with regular waves, is not great. I joined the figures, and my matlab code below. Maybe you will be able to see what is wrong.

Wsum=zeros(50001,1);
Hsum=zeros(50001,1);
for i=1:1:6;
nb=num2str(i)
name = [Directory ‘\white_noise_’ nb]
fid = fopen(name, ‘rt’);
for J=1:1:8
fgetl(fid);
end
M1 = fscanf(fid, ‘%f’, [10,inf]).’;
Wave1Elev1=M1(:,4) % get the Wave elevation time series
WRPHeave1=M1(:,3) % get the heave displacement time series
fclose(fid);

Wsum = Wsum + Wave1Elev1/6
Hsum = Hsum + WRPHeave1/6
end

Fs=1/0.05 % in time series, timestep= 0.05s
N = length(Wsum);
F = fft(Wsum);
F = F(1:N/2+1);
ASD = (1/(FsN)) * abs(F).^2;
ASD(2:end-1) = 2
ASD(2:end-1);
freq = 0:Fs/length(Wsum):Fs/2;
subplot(311)
plot(freq, ASD) %plot the auto-spectral density of Wsum
legend(‘auto-spectral density of Wsum’)

F1 = fft(Wsum);
F2 = fft(Hsum);
F1 = F1(1:N/2+1);
F2 = F2(1:N/2+1);
CSD = (1/(FsN))abs( F2 . conj(F1));
CSD(2:end-1) = 2
CSD(2:end-1)
subplot(312)
plot(freq, CSD) %plot the cross spectral density of Wsum and Hsum
legend('cross spectral density of Wsum and Hsum ')

RAO = CSD./ASD
subplot(313)
plot(freq, RAO) %plot the RAO
legend(‘RAO’)

Best regards,

Julie Marchebois

Dear Julie,

I haven’t tried to review your code to understand why your results are not what you’d expect, but a MATLAB script (RAO.m) for calculating RAOs based on the approach outlined in our ISOPE paper is available in my post dated Mar 09, 2018 in the following forum topic: http://forums.nrel.gov/t/rao-oc3hywind-turbine/1753/1.

Best regards,

Dear Mr,

Thank you for the advice. I used your matlab code RAO.m and here is what I get.

There is still a problem for low frequencies. I would like to know if it could come from my .hst file. Indeed, I didn’t use wamit to obtain it and I calculated the hydrostatic coefficients myself by using the definitions of the wamit documentation : https://www.wamit.com/manual6.4/Chap4.pdf.

Best regards,

Julie Marchebois

Dear Julie,

Can you clarify what your concern is? Are you concerned with the spike in PSD at zero frequency? This results from some nonzero mean heave displacement.

Best regards,

Dear Jason,

Thank you for your answers, I finally get a RAO of heave motion that seems reasonable to me with a RAO, which tends to one in low frequencies. Now I would like to try to output an RAO of heave motion with a Jonswap waves spectrum, by using the wavemod 2. I read that it wasn’t easy to get this. Anyway I launched a simulation with an Hs of 3.5m and a Tp of 15s. Then I computed the cross spectral density with the help of cspd, a matlab function, and the auto spectral density with the matlab function pwelch. Then, I plot the cross spectal density divided by the auto-spectral density and here is what I get.
I don’t know what to think of this and I would like to know if there is a way to improve my result. I saw in several topics that you advise people to compute the psd with the code “PSD.m”, but I guess we can’t use this code to calculate the CPSD too ?

Best regards,
Julie Marchebois

Dear Julio,

I’m not familiar with these MATLAB functions to know if they do something similar, but the Jason_PSD.m script I’ve provided computes the PSD directly via FFT and the RAO.m script I’ve shared also computes the auto-spectral density and cross-spectral density via FFT. You can review the scripts to see the details.

The big problem with using a JONSWAP spectrum to compute RAOs is that the the JONSWAP spectrum has little energy outside the primary range of wave frequencies and the natural frequencies of the floater likely fall outside of this frequency range.

Best regards,