Skip to content

Commit 506785b

Browse files
jtgrasbjtgrasb
committed
Add full directional wave examples and test
1 parent a8d0434 commit 506785b

11 files changed

Lines changed: 248 additions & 0 deletions

‎Full_Directional_Waves/MHKiT_NDBC_Waves/OSWEC.slx‎

26.5 KB
Binary file not shown.

‎Full_Directional_Waves/MHKiT_NDBC_Waves/fullDirSpectrum.mat‎

63.2 KB
Binary file not shown.

‎Full_Directional_Waves/MHKiT_NDBC_Waves/userDefinedFunctions.m‎

Lines changed: 22 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,22 @@
1+
%Example of user input MATLAB file for post processing
2+
3+
%Plot waves
4+
waves.plotElevation(simu.rampTime);
5+
try
6+
waves.plotSpectrum();
7+
catch
8+
end
9+
10+
% Plot RY forces for body 1
11+
output.plotForces(1,5)
12+
13+
%Plot RY response for body 1
14+
output.plotResponse(1,5);
15+
16+
% Plot x forces for body 2
17+
output.plotForces(2,1)
18+
19+
% Save waves and response as video
20+
% output.saveViz(simu,body,waves,...
21+
% 'timesPerFrame',5,'axisLimits',[-50 50 -50 50 -12 20],...
22+
% 'startEndTime',[100 125]);

‎Full_Directional_Waves/MHKiT_NDBC_Waves/wecSimInputFile.m‎

Lines changed: 41 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,41 @@
1+
%% Simulation Data
2+
simu = simulationClass(); % Initialize Simulation Class
3+
simu.simMechanicsFile = 'OSWEC.slx'; % Specify Simulink Model File
4+
simu.mode = 'normal'; % Specify Simulation Mode ('normal','accelerator','rapid-accelerator')
5+
simu.explorer = 'on'; % Turn SimMechanics Explorer (on/off)
6+
simu.startTime = 0; % Simulation Start Time [s]
7+
simu.endTime = 400; % Simulation End Time [s]
8+
simu.rampTime = 100; % Wave Ramp Time [s]
9+
simu.solver = 'ode4'; % simu.solver = 'ode4' for fixed step & simu.solver = 'ode45' for variable step
10+
simu.dt = 0.05; % Simulation time-step [s]
11+
simu.cicEndTime = 30;
12+
13+
%% Wave Information
14+
% Full directional waves
15+
waves = waveClass('spectrumImportFullDir');
16+
waves.spectrumFile = ('fullDirSpectrum.mat');
17+
waves.phaseSeed = 1;
18+
19+
%% Body Data
20+
% Flap
21+
body(1) = bodyClass('../../_Common_Input_Files/OSWEC/hydroData/oswec.h5'); % Initialize bodyClass for Flap
22+
body(1).geometryFile = '../../_Common_Input_Files/OSWEC/geometry/flap.stl'; % Geometry File
23+
body(1).mass = 127000; % User-Defined mass [kg]
24+
body(1).inertia = [1.85e6 1.85e6 1.85e6]; % Moment of Inertia [kg-m^2]
25+
26+
% Base
27+
body(2) = bodyClass('../../_Common_Input_Files/OSWEC/hydroData/oswec.h5'); % Initialize bodyClass for Base
28+
body(2).geometryFile = '../../_Common_Input_Files/OSWEC/geometry/base.stl'; % Geometry File
29+
body(2).mass = 999; % Placeholder mass for a fixed body
30+
body(2).inertia = [999 999 999]; % Placeholder inertia for a fixed body
31+
32+
%% PTO and Constraint Parameters
33+
% Fixed
34+
constraint(1)= constraintClass('Constraint1'); % Initialize ConstraintClass for Constraint1
35+
constraint(1).location = [0 0 -10]; % Constraint Location [m]
36+
37+
% Rotational PTO
38+
pto(1) = ptoClass('PTO1'); % Initialize ptoClass for PTO1
39+
pto(1).stiffness = 0; % PTO Stiffness Coeff [Nm/rad]
40+
pto(1).damping = 12000; % PTO Damping Coeff [Nsm/rad]
41+
pto(1).location = [0 0 -8.9]; % PTO Location [m]

‎Full_Directional_Waves/OOI_Waves/OSWEC.slx‎

26.5 KB
Binary file not shown.

‎Full_Directional_Waves/OOI_Waves/convertOOIToIEC.m‎

Lines changed: 51 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,51 @@
1+
function [frequenciesIEC, spectrumIEC, spreadIEC, directionsIEC] = convertOOIToIEC(spectrumDataOOI, directionBins, plotFlag)
2+
3+
directionBins = wrapTo180(directionBins);
4+
if any(directionBins == 180) && any(directionBins == -180)
5+
warning('Directions include both -180 and 180. Removing -180 to avoid duplicate directions.')
6+
directionBins(directionBins == -180) = [];
7+
end
8+
if length(unique(mod(directionBins + 180, 360) - 180)) < length(directionBins)
9+
error('Duplicate directions were found.')
10+
end
11+
12+
frequenciesOOI = spectrumDataOOI(:,1);
13+
spectrumOOI = spectrumDataOOI(:,2);
14+
directionsOOI = spectrumDataOOI(:,3);
15+
spreadOOI = spectrumDataOOI(:,4);
16+
17+
% convert to IEC:
18+
directionsIEC = directionBins(:);
19+
frequenciesIEC = frequenciesOOI;
20+
spectrumIEC = spectrumOOI;
21+
spreadIEC = zeros(length(frequenciesIEC),length(directionsIEC));
22+
23+
% for spread, we need to apply gaussian distribution
24+
for ii = 1:length(frequenciesOOI)
25+
% convert to IEC spread calculation
26+
directions = deg2rad(wrapTo180(directionBins - wrapTo180(directionsOOI(ii))));
27+
energySpread = (1./(deg2rad(spreadOOI(ii)).*sqrt(2*pi))) .* exp(-(directions.^2) ./ (2.*deg2rad(spreadOOI(ii)).^2));
28+
checkSum = trapz(deg2rad(directionBins),energySpread);
29+
if checkSum < 0.95 % if this is true, then less than 95% of the initial energy at this frequency is contained over the considered directions.
30+
warning('Number of directional bins inadequate at frequency number %i. Directional approximation weak. \n \r', ii);
31+
end
32+
energySpread = energySpread./checkSum;
33+
spreadIEC(ii,:) = energySpread;
34+
end
35+
36+
if plotFlag == 1
37+
spectrumFullDir = spreadIEC.*spectrumIEC;
38+
meanDirection = sum(wrapTo180(directionsOOI).*spectrumOOI)/sum(spectrumOOI);
39+
40+
[plotDirs,plotFreqs] = meshgrid(directionsIEC,frequenciesIEC);
41+
42+
figure()
43+
polarscatter(deg2rad(plotDirs(:)),plotFreqs(:),4,spectrumFullDir(:),'filled')
44+
hold on
45+
polarplot(deg2rad([meanDirection meanDirection]), [0 max(plotFreqs(:))], 'k--'); %
46+
c = colorbar;
47+
c.Label.String = 'Spectrum (m^2/Hz/deg)';
48+
title('Elevation variance spectrum');
49+
legend('spectrum','approx. mean direction','Location','northwest')
50+
end
51+
end

‎Full_Directional_Waves/OOI_Waves/dirSpectrumOOI.mat‎

1.18 KB
Binary file not shown.

‎Full_Directional_Waves/OOI_Waves/fullDirSpectrum.mat‎

88.3 KB
Binary file not shown.

‎Full_Directional_Waves/OOI_Waves/userDefinedFunctions.m‎

Lines changed: 22 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,22 @@
1+
%Example of user input MATLAB file for post processing
2+
3+
%Plot waves
4+
waves.plotElevation(simu.rampTime);
5+
try
6+
waves.plotSpectrum();
7+
catch
8+
end
9+
10+
% Plot RY forces for body 1
11+
output.plotForces(1,5)
12+
13+
%Plot RY response for body 1
14+
output.plotResponse(1,5);
15+
16+
% Plot x forces for body 2
17+
output.plotForces(2,1)
18+
19+
% Save waves and response as video
20+
% output.saveViz(simu,body,waves,...
21+
% 'timesPerFrame',5,'axisLimits',[-50 50 -50 50 -12 20],...
22+
% 'startEndTime',[100 125]);

‎Full_Directional_Waves/OOI_Waves/wecSimInputFile.m‎

Lines changed: 48 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,48 @@
1+
%% Simulation Data
2+
simu = simulationClass(); % Initialize Simulation Class
3+
simu.simMechanicsFile = 'OSWEC.slx'; % Specify Simulink Model File
4+
simu.mode = 'normal'; % Specify Simulation Mode ('normal','accelerator','rapid-accelerator')
5+
simu.explorer = 'on'; % Turn SimMechanics Explorer (on/off)
6+
simu.startTime = 0; % Simulation Start Time [s]
7+
simu.endTime = 400; % Simulation End Time [s]
8+
simu.rampTime = 100; % Wave Ramp Time [s]
9+
simu.solver = 'ode4'; % simu.solver = 'ode4' for fixed step & simu.solver = 'ode45' for variable step
10+
simu.dt = 0.05; % Simulation time-step [s]
11+
simu.cicEndTime = 30;
12+
13+
%% Wave Information
14+
load dirSpectrumOOI.mat
15+
spectrumDataOOI = dataWaveSnip;
16+
directions = -179:2:179;
17+
18+
[frequencies, spectrum, spread, directions] = convertOOIToIEC(spectrumDataOOI, directions, 1);
19+
20+
save 'fullDirSpectrum.mat' spectrum spread frequencies directions
21+
22+
waves = waveClass('spectrumImportFullDir');
23+
waves.spectrumFile = ('fullDirSpectrum.mat');
24+
waves.phaseSeed = 1;
25+
26+
%% Body Data
27+
% Flap
28+
body(1) = bodyClass('../../_Common_Input_Files/OSWEC/hydroData/oswec.h5'); % Initialize bodyClass for Flap
29+
body(1).geometryFile = '../../_Common_Input_Files/OSWEC/geometry/flap.stl'; % Geometry File
30+
body(1).mass = 127000; % User-Defined mass [kg]
31+
body(1).inertia = [1.85e6 1.85e6 1.85e6]; % Moment of Inertia [kg-m^2]
32+
33+
% Base
34+
body(2) = bodyClass('../../_Common_Input_Files/OSWEC/hydroData/oswec.h5'); % Initialize bodyClass for Base
35+
body(2).geometryFile = '../../_Common_Input_Files/OSWEC/geometry/base.stl'; % Geometry File
36+
body(2).mass = 999; % Placeholder mass for a fixed body
37+
body(2).inertia = [999 999 999]; % Placeholder inertia for a fixed body
38+
39+
%% PTO and Constraint Parameters
40+
% Fixed
41+
constraint(1)= constraintClass('Constraint1'); % Initialize ConstraintClass for Constraint1
42+
constraint(1).location = [0 0 -10]; % Constraint Location [m]
43+
44+
% Rotational PTO
45+
pto(1) = ptoClass('PTO1'); % Initialize ptoClass for PTO1
46+
pto(1).stiffness = 0; % PTO Stiffness Coeff [Nm/rad]
47+
pto(1).damping = 12000; % PTO Damping Coeff [Nsm/rad]
48+
pto(1).location = [0 0 -8.9]; % PTO Location [m]

0 commit comments

Comments
 (0)