Skip to content

Commit 9a7a98c

Browse files
jtgrasbjtgrasb
authored
Add full directional wave examples and test (WEC-Sim#103)
* Add full directional wave examples and test * Combine full directional wave applications * remove created file * use R2024b for changed tests * Add products.txt
1 parent b7e7ab4 commit 9a7a98c

9 files changed

Lines changed: 189 additions & 1 deletion

File tree

‎.github/workflows/changed-tests.yml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -37,7 +37,7 @@ jobs:
3737
- name: Install MATLAB
3838
uses: matlab-actions/setup-matlab@v2
3939
with:
40-
release: R2024a
40+
release: R2024b
4141
- name: Get test target folder
4242
uses: matlab-actions/run-command@v2
4343
with:

‎Full_Directional_Waves/OSWEC.slx‎

26.5 KB
Binary file not shown.

‎Full_Directional_Waves/TestFullDirectionalWaves.m‎

Lines changed: 57 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,57 @@
1+
classdef TestFullDirectionalWaves < matlab.unittest.TestCase
2+
3+
properties
4+
OriginalDefault
5+
testDir
6+
h5Dir = '../_Common_Input_Files/OSWEC/hydroData/'
7+
h5Name = 'oswec.h5'
8+
outName = 'oswec.out'
9+
end
10+
11+
methods (Access = 'public')
12+
function obj = TestFullDirectionalWaves
13+
obj.testDir = fileparts(mfilename('fullpath'));
14+
end
15+
end
16+
17+
methods (TestMethodSetup)
18+
function killPlots (~)
19+
set(0,'DefaultFigureVisible','off');
20+
end
21+
end
22+
23+
methods(TestClassSetup)
24+
function captureVisibility(testCase)
25+
testCase.OriginalDefault = get(0,'DefaultFigureVisible');
26+
end
27+
function runBemio(testCase)
28+
cd(testCase.h5Dir);
29+
if isfile(testCase.h5Name)
30+
fprintf('runBemio skipped, *.h5 already exists\n')
31+
else
32+
bemio
33+
end
34+
cd(testCase.testDir)
35+
end
36+
end
37+
38+
methods(TestMethodTeardown)
39+
function returnHome(testCase)
40+
cd(testCase.testDir)
41+
end
42+
end
43+
44+
methods(TestClassTeardown)
45+
function checkVisibilityRestored(testCase)
46+
set(0,'DefaultFigureVisible',testCase.OriginalDefault);
47+
testCase.assertEqual(get(0,'DefaultFigureVisible'), ...
48+
testCase.OriginalDefault);
49+
end
50+
end
51+
52+
methods(Test)
53+
function testFullDirectionalWaves(testCase)
54+
wecSim
55+
end
56+
end
57+
end

‎Full_Directional_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/dirSpectrumOOI.mat‎

1.18 KB
Binary file not shown.

‎Full_Directional_Waves/fullDirSpectrum.mat‎

63.2 KB
Binary file not shown.

‎Full_Directional_Waves/products.txt‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,3 @@
1+
Simulink
2+
Simscape
3+
Simscape_Multibody

‎Full_Directional_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/wecSimInputFile.m‎

Lines changed: 55 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,55 @@
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+
15+
% Full directional waves base case
16+
waves = waveClass('spectrumImportFullDir');
17+
waves.spectrumFile = ('fullDirSpectrum.mat');
18+
waves.phaseSeed = 1;
19+
20+
% % Uncomment this section to run full directional waves with conversion from OOI
21+
% load dirSpectrumOOI.mat
22+
% spectrumDataOOI = dataWaveSnip;
23+
% directions = -179:2:179;
24+
%
25+
% [frequencies, spectrum, spread, directions] = convertOOIToIEC(spectrumDataOOI, directions, 1);
26+
%
27+
% save 'fullDirSpectrumOOI.mat' spectrum spread frequencies directions
28+
%
29+
% waves = waveClass('spectrumImportFullDir');
30+
% waves.spectrumFile = ('fullDirSpectrumOOI.mat');
31+
% waves.phaseSeed = 1;
32+
33+
%% Body Data
34+
% Flap
35+
body(1) = bodyClass('../_Common_Input_Files/OSWEC/hydroData/oswec.h5'); % Initialize bodyClass for Flap
36+
body(1).geometryFile = '../_Common_Input_Files/OSWEC/geometry/flap.stl'; % Geometry File
37+
body(1).mass = 127000; % User-Defined mass [kg]
38+
body(1).inertia = [1.85e6 1.85e6 1.85e6]; % Moment of Inertia [kg-m^2]
39+
40+
% Base
41+
body(2) = bodyClass('../_Common_Input_Files/OSWEC/hydroData/oswec.h5'); % Initialize bodyClass for Base
42+
body(2).geometryFile = '../_Common_Input_Files/OSWEC/geometry/base.stl'; % Geometry File
43+
body(2).mass = 999; % Placeholder mass for a fixed body
44+
body(2).inertia = [999 999 999]; % Placeholder inertia for a fixed body
45+
46+
%% PTO and Constraint Parameters
47+
% Fixed
48+
constraint(1)= constraintClass('Constraint1'); % Initialize ConstraintClass for Constraint1
49+
constraint(1).location = [0 0 -10]; % Constraint Location [m]
50+
51+
% Rotational PTO
52+
pto(1) = ptoClass('PTO1'); % Initialize ptoClass for PTO1
53+
pto(1).stiffness = 0; % PTO Stiffness Coeff [Nm/rad]
54+
pto(1).damping = 12000; % PTO Damping Coeff [Nsm/rad]
55+
pto(1).location = [0 0 -8.9]; % PTO Location [m]

0 commit comments

Comments
 (0)