diff --git a/Full_Directional_Waves/MHKiT_NDBC_Waves/OSWEC.slx b/Full_Directional_Waves/MHKiT_NDBC_Waves/OSWEC.slx new file mode 100644 index 00000000..42d0e230 Binary files /dev/null and b/Full_Directional_Waves/MHKiT_NDBC_Waves/OSWEC.slx differ diff --git a/Full_Directional_Waves/MHKiT_NDBC_Waves/fullDirSpectrum.mat b/Full_Directional_Waves/MHKiT_NDBC_Waves/fullDirSpectrum.mat new file mode 100644 index 00000000..315cff9d Binary files /dev/null and b/Full_Directional_Waves/MHKiT_NDBC_Waves/fullDirSpectrum.mat differ diff --git a/Full_Directional_Waves/MHKiT_NDBC_Waves/userDefinedFunctions.m b/Full_Directional_Waves/MHKiT_NDBC_Waves/userDefinedFunctions.m new file mode 100644 index 00000000..cbaeffcf --- /dev/null +++ b/Full_Directional_Waves/MHKiT_NDBC_Waves/userDefinedFunctions.m @@ -0,0 +1,22 @@ +%Example of user input MATLAB file for post processing + +%Plot waves +waves.plotElevation(simu.rampTime); +try + waves.plotSpectrum(); +catch +end + +% Plot RY forces for body 1 +output.plotForces(1,5) + +%Plot RY response for body 1 +output.plotResponse(1,5); + +% Plot x forces for body 2 +output.plotForces(2,1) + +% Save waves and response as video +% output.saveViz(simu,body,waves,... +% 'timesPerFrame',5,'axisLimits',[-50 50 -50 50 -12 20],... +% 'startEndTime',[100 125]); diff --git a/Full_Directional_Waves/MHKiT_NDBC_Waves/wecSimInputFile.m b/Full_Directional_Waves/MHKiT_NDBC_Waves/wecSimInputFile.m new file mode 100644 index 00000000..9271ae59 --- /dev/null +++ b/Full_Directional_Waves/MHKiT_NDBC_Waves/wecSimInputFile.m @@ -0,0 +1,41 @@ +%% Simulation Data +simu = simulationClass(); % Initialize Simulation Class +simu.simMechanicsFile = 'OSWEC.slx'; % Specify Simulink Model File +simu.mode = 'normal'; % Specify Simulation Mode ('normal','accelerator','rapid-accelerator') +simu.explorer = 'on'; % Turn SimMechanics Explorer (on/off) +simu.startTime = 0; % Simulation Start Time [s] +simu.endTime = 400; % Simulation End Time [s] +simu.rampTime = 100; % Wave Ramp Time [s] +simu.solver = 'ode4'; % simu.solver = 'ode4' for fixed step & simu.solver = 'ode45' for variable step +simu.dt = 0.05; % Simulation time-step [s] +simu.cicEndTime = 30; + +%% Wave Information +% Full directional waves +waves = waveClass('spectrumImportFullDir'); +waves.spectrumFile = ('fullDirSpectrum.mat'); +waves.phaseSeed = 1; + +%% Body Data +% Flap +body(1) = bodyClass('../../_Common_Input_Files/OSWEC/hydroData/oswec.h5'); % Initialize bodyClass for Flap +body(1).geometryFile = '../../_Common_Input_Files/OSWEC/geometry/flap.stl'; % Geometry File +body(1).mass = 127000; % User-Defined mass [kg] +body(1).inertia = [1.85e6 1.85e6 1.85e6]; % Moment of Inertia [kg-m^2] + +% Base +body(2) = bodyClass('../../_Common_Input_Files/OSWEC/hydroData/oswec.h5'); % Initialize bodyClass for Base +body(2).geometryFile = '../../_Common_Input_Files/OSWEC/geometry/base.stl'; % Geometry File +body(2).mass = 999; % Placeholder mass for a fixed body +body(2).inertia = [999 999 999]; % Placeholder inertia for a fixed body + +%% PTO and Constraint Parameters +% Fixed +constraint(1)= constraintClass('Constraint1'); % Initialize ConstraintClass for Constraint1 +constraint(1).location = [0 0 -10]; % Constraint Location [m] + +% Rotational PTO +pto(1) = ptoClass('PTO1'); % Initialize ptoClass for PTO1 +pto(1).stiffness = 0; % PTO Stiffness Coeff [Nm/rad] +pto(1).damping = 12000; % PTO Damping Coeff [Nsm/rad] +pto(1).location = [0 0 -8.9]; % PTO Location [m] diff --git a/Full_Directional_Waves/OOI_Waves/OSWEC.slx b/Full_Directional_Waves/OOI_Waves/OSWEC.slx new file mode 100644 index 00000000..42d0e230 Binary files /dev/null and b/Full_Directional_Waves/OOI_Waves/OSWEC.slx differ diff --git a/Full_Directional_Waves/OOI_Waves/convertOOIToIEC.m b/Full_Directional_Waves/OOI_Waves/convertOOIToIEC.m new file mode 100644 index 00000000..0f19c4c4 --- /dev/null +++ b/Full_Directional_Waves/OOI_Waves/convertOOIToIEC.m @@ -0,0 +1,51 @@ +function [frequenciesIEC, spectrumIEC, spreadIEC, directionsIEC] = convertOOIToIEC(spectrumDataOOI, directionBins, plotFlag) + + directionBins = wrapTo180(directionBins); + if any(directionBins == 180) && any(directionBins == -180) + warning('Directions include both -180 and 180. Removing -180 to avoid duplicate directions.') + directionBins(directionBins == -180) = []; + end + if length(unique(mod(directionBins + 180, 360) - 180)) < length(directionBins) + error('Duplicate directions were found.') + end + + frequenciesOOI = spectrumDataOOI(:,1); + spectrumOOI = spectrumDataOOI(:,2); + directionsOOI = spectrumDataOOI(:,3); + spreadOOI = spectrumDataOOI(:,4); + + % convert to IEC: + directionsIEC = directionBins(:); + frequenciesIEC = frequenciesOOI; + spectrumIEC = spectrumOOI; + spreadIEC = zeros(length(frequenciesIEC),length(directionsIEC)); + + % for spread, we need to apply gaussian distribution + for ii = 1:length(frequenciesOOI) + % convert to IEC spread calculation + directions = deg2rad(wrapTo180(directionBins - wrapTo180(directionsOOI(ii)))); + energySpread = (1./(deg2rad(spreadOOI(ii)).*sqrt(2*pi))) .* exp(-(directions.^2) ./ (2.*deg2rad(spreadOOI(ii)).^2)); + checkSum = trapz(deg2rad(directionBins),energySpread); + 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. + warning('Number of directional bins inadequate at frequency number %i. Directional approximation weak. \n \r', ii); + end + energySpread = energySpread./checkSum; + spreadIEC(ii,:) = energySpread; + end + + if plotFlag == 1 + spectrumFullDir = spreadIEC.*spectrumIEC; + meanDirection = sum(wrapTo180(directionsOOI).*spectrumOOI)/sum(spectrumOOI); + + [plotDirs,plotFreqs] = meshgrid(directionsIEC,frequenciesIEC); + + figure() + polarscatter(deg2rad(plotDirs(:)),plotFreqs(:),4,spectrumFullDir(:),'filled') + hold on + polarplot(deg2rad([meanDirection meanDirection]), [0 max(plotFreqs(:))], 'k--'); % + c = colorbar; + c.Label.String = 'Spectrum (m^2/Hz/deg)'; + title('Elevation variance spectrum'); + legend('spectrum','approx. mean direction','Location','northwest') + end +end \ No newline at end of file diff --git a/Full_Directional_Waves/OOI_Waves/dirSpectrumOOI.mat b/Full_Directional_Waves/OOI_Waves/dirSpectrumOOI.mat new file mode 100644 index 00000000..85f91567 Binary files /dev/null and b/Full_Directional_Waves/OOI_Waves/dirSpectrumOOI.mat differ diff --git a/Full_Directional_Waves/OOI_Waves/fullDirSpectrum.mat b/Full_Directional_Waves/OOI_Waves/fullDirSpectrum.mat new file mode 100644 index 00000000..498f6b54 Binary files /dev/null and b/Full_Directional_Waves/OOI_Waves/fullDirSpectrum.mat differ diff --git a/Full_Directional_Waves/OOI_Waves/userDefinedFunctions.m b/Full_Directional_Waves/OOI_Waves/userDefinedFunctions.m new file mode 100644 index 00000000..cbaeffcf --- /dev/null +++ b/Full_Directional_Waves/OOI_Waves/userDefinedFunctions.m @@ -0,0 +1,22 @@ +%Example of user input MATLAB file for post processing + +%Plot waves +waves.plotElevation(simu.rampTime); +try + waves.plotSpectrum(); +catch +end + +% Plot RY forces for body 1 +output.plotForces(1,5) + +%Plot RY response for body 1 +output.plotResponse(1,5); + +% Plot x forces for body 2 +output.plotForces(2,1) + +% Save waves and response as video +% output.saveViz(simu,body,waves,... +% 'timesPerFrame',5,'axisLimits',[-50 50 -50 50 -12 20],... +% 'startEndTime',[100 125]); diff --git a/Full_Directional_Waves/OOI_Waves/wecSimInputFile.m b/Full_Directional_Waves/OOI_Waves/wecSimInputFile.m new file mode 100644 index 00000000..90bef2dc --- /dev/null +++ b/Full_Directional_Waves/OOI_Waves/wecSimInputFile.m @@ -0,0 +1,48 @@ +%% Simulation Data +simu = simulationClass(); % Initialize Simulation Class +simu.simMechanicsFile = 'OSWEC.slx'; % Specify Simulink Model File +simu.mode = 'normal'; % Specify Simulation Mode ('normal','accelerator','rapid-accelerator') +simu.explorer = 'on'; % Turn SimMechanics Explorer (on/off) +simu.startTime = 0; % Simulation Start Time [s] +simu.endTime = 400; % Simulation End Time [s] +simu.rampTime = 100; % Wave Ramp Time [s] +simu.solver = 'ode4'; % simu.solver = 'ode4' for fixed step & simu.solver = 'ode45' for variable step +simu.dt = 0.05; % Simulation time-step [s] +simu.cicEndTime = 30; + +%% Wave Information +load dirSpectrumOOI.mat +spectrumDataOOI = dataWaveSnip; +directions = -179:2:179; + +[frequencies, spectrum, spread, directions] = convertOOIToIEC(spectrumDataOOI, directions, 1); + +save 'fullDirSpectrum.mat' spectrum spread frequencies directions + +waves = waveClass('spectrumImportFullDir'); +waves.spectrumFile = ('fullDirSpectrum.mat'); +waves.phaseSeed = 1; + +%% Body Data +% Flap +body(1) = bodyClass('../../_Common_Input_Files/OSWEC/hydroData/oswec.h5'); % Initialize bodyClass for Flap +body(1).geometryFile = '../../_Common_Input_Files/OSWEC/geometry/flap.stl'; % Geometry File +body(1).mass = 127000; % User-Defined mass [kg] +body(1).inertia = [1.85e6 1.85e6 1.85e6]; % Moment of Inertia [kg-m^2] + +% Base +body(2) = bodyClass('../../_Common_Input_Files/OSWEC/hydroData/oswec.h5'); % Initialize bodyClass for Base +body(2).geometryFile = '../../_Common_Input_Files/OSWEC/geometry/base.stl'; % Geometry File +body(2).mass = 999; % Placeholder mass for a fixed body +body(2).inertia = [999 999 999]; % Placeholder inertia for a fixed body + +%% PTO and Constraint Parameters +% Fixed +constraint(1)= constraintClass('Constraint1'); % Initialize ConstraintClass for Constraint1 +constraint(1).location = [0 0 -10]; % Constraint Location [m] + +% Rotational PTO +pto(1) = ptoClass('PTO1'); % Initialize ptoClass for PTO1 +pto(1).stiffness = 0; % PTO Stiffness Coeff [Nm/rad] +pto(1).damping = 12000; % PTO Damping Coeff [Nsm/rad] +pto(1).location = [0 0 -8.9]; % PTO Location [m] diff --git a/Full_Directional_Waves/TestFullDirectionalWaves.m b/Full_Directional_Waves/TestFullDirectionalWaves.m new file mode 100644 index 00000000..65534c8b --- /dev/null +++ b/Full_Directional_Waves/TestFullDirectionalWaves.m @@ -0,0 +1,64 @@ +classdef TestFullDirectionalWaves < matlab.unittest.TestCase + + properties + OriginalDefault + testDir + h5Dir = '../_Common_Input_Files/OSWEC/hydroData/' + h5Name = 'oswec.h5' + outName = 'oswec.out' + end + + methods (Access = 'public') + function obj = TestFullDirectionalWaves + obj.testDir = fileparts(mfilename('fullpath')); + end + end + + methods (TestMethodSetup) + function killPlots (~) + set(0,'DefaultFigureVisible','off'); + end + end + + methods(TestClassSetup) + function captureVisibility(testCase) + testCase.OriginalDefault = get(0,'DefaultFigureVisible'); + end + function runBemio(testCase) + cd(testCase.h5Dir); + if isfile(testCase.h5Name) + fprintf('runBemio skipped, *.h5 already exists\n') + else + bemio + end + cd(testCase.testDir) + end + end + + methods(TestMethodTeardown) + function returnHome(testCase) + cd(testCase.testDir) + end + end + + methods(TestClassTeardown) + function checkVisibilityRestored(testCase) + set(0,'DefaultFigureVisible',testCase.OriginalDefault); + testCase.assertEqual(get(0,'DefaultFigureVisible'), ... + testCase.OriginalDefault); + end + end + + methods(Test) + function testMHKiT_NDBC_Waves(testCase) + cd MHKiT_NDBC_Waves + wecSim + close_system('OSWEC',0) + end + function testOOI_Waves(testCase) + cd OOI_Waves + wecSim + close_system('OSWEC',0) + end + end +end