Skip to content

Commit 34e61e7

Browse files
akeesteH0R5E
akeeste
and
H0R5E
authored
Variable hydro demonstrations (#65)
* initial test cast - variable hydro vs passive yaw * initial passive yaw comparison * working case and scripted comparison to passive yaw simulations * clean up bemio function * increase BEM resolution * increase BEM resolution for irregular case * add variable hydro test * fix typo in readme * Move products file to top level directory * default to and test regular wave cases for speed purposes --------- Co-authored-by: Mathew Topper <dataonlygreater@gmail.com>
1 parent 8c66cfa commit 34e61e7

13 files changed

Lines changed: 362 additions & 1 deletion

File tree

‎Cable/README.md‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -10,4 +10,4 @@
1010

1111
**Note:** More WEC information can be found: https://link.springer.com/article/10.1007/s40722-021-00197-9
1212

13-
Example using WEC-Sim to simulate a [Cable](http://wec-sim.github.io/WEC-Sim/advanced_features.html) connecting two rigid bodies for the ??? geometry.
13+
Example using WEC-Sim to simulate a [Cable](http://wec-sim.github.io/WEC-Sim/advanced_features.html) connecting two rigid bodies for the MBARI geometry.

‎Variable_Hydro/Passive_Yaw/OSWEC.slx‎

41.4 KB
Binary file not shown.

‎Variable_Hydro/Passive_Yaw/README.md‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,11 @@
1+
# Cable
2+
3+
**Author:** Adam Keester
4+
5+
**Version:** WEC-Sim v6.1
6+
7+
**Geometry:** OSWEC
8+
9+
**Dependency:** N/A
10+
11+
Example demonstrating WEC-Sim's variable hydrodynamics feature by comparing it to the passive yaw implementation.

‎Variable_Hydro/Passive_Yaw/TestVariableHydro.m‎

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

‎Variable_Hydro/Passive_Yaw/calcIndex.m‎

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,30 @@
1+
function body1_hydroForceIndex = calcIndex(Rz, waveDirection, bemDirections, previousIndex)
2+
% This case is mimicking passive yaw. Passive yaw is going to interpolate
3+
% and pull BEM data for the incident wave direction closest to the relative
4+
% angle between the incoming wave and the device yaw position. So with an
5+
% incoming wave angle of 10 deg, and some instantaneous yaw position:
6+
%
7+
% yaw = -10 --> BEM at 20 deg
8+
% yaw = 0 --> BEM at 10 deg
9+
% yaw = 10 --> BEM at 0 deg
10+
% yaw = 20 --> BEM at -10 deg
11+
12+
relativeAngle = waveDirection - Rz*180/pi;
13+
14+
% % This indexing method is indentical to the passive yaw implementation
15+
% % but converges more slowly. To use, add a `memory` block in the simulink
16+
% % model that connects body1_hydroForceIndex to previousIndex.
17+
% previousDir = bemDirections(previousIndex);
18+
% dTheta = bemDirections(2) - bemDirections(1); % assume BEM discretization is constant and equivalent to passive yaw's threshold
19+
% if abs(relativeAngle - previousDir) > dTheta
20+
% [~, body1_hydroForceIndex] = min(abs(bemDirections-relativeAngle));
21+
% else
22+
% body1_hydroForceIndex = previousIndex;
23+
% end
24+
25+
% This indexing method (always choosing the closest BEM dataset) is not
26+
% identical to the passive yaw implementation, but is simple, precise, and
27+
% converges quickly.
28+
[~, body1_hydroForceIndex] = min(abs(bemDirections-relativeAngle));
29+
30+
end

‎Variable_Hydro/Passive_Yaw/hydroData/bemio.m‎

Lines changed: 78 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,78 @@
1+
%% Read OSWEC hydro data
2+
hydro = struct();
3+
hydro = readWAMIT(hydro,'../../../_Common_Input_Files/OSWEC/hydroData/oswec.out',[]);
4+
hydro = radiationIRF(hydro,40,[],[],[],[]);
5+
hydro = excitationIRF(hydro,40,[],[],[],[]);
6+
hydro.plotDirections = [36 1 2 3];
7+
hydro.plotDofs = [6 6];
8+
hydro.plotBodies = 1;
9+
% plotBEMIO(hydro)
10+
% hydro = writeBEMIOH5(hydro);
11+
12+
%% Split hydro into multiple hydro structures by wave direction
13+
iDirs = [32 33 34 35 36 1 2 3 4 5 6]; % select direction indices to use
14+
for i = 1:length(iDirs)
15+
iDir = iDirs(i);
16+
17+
temp_hydro = hydro;
18+
vars = {'ex_K','ex_ma', 'ex_ph', 'ex_re', 'ex_im', ...
19+
'sc_ma', 'sc_ph', 'sc_re', 'sc_im', ...
20+
'fk_ma', 'fk_ph', 'fk_re', 'fk_im'}; % directionally dependent variables
21+
22+
for iVar = 1:length(vars)
23+
temp_hydro.(vars{iVar}) = hydro.(vars{iVar})(:,iDir,:);
24+
end
25+
temp_hydro.Nh = 1;
26+
% temp_hydro.file = [hydro.file '_' num2str(temp_hydro.theta)];
27+
temp_hydro.file = [hydro.file '_' num2str(hydro.theta(iDir))];
28+
temp_hydro.theta = 10;
29+
30+
% theta(i) = temp_hydro.theta;
31+
theta(i) = hydro.theta(iDir);
32+
33+
% Assign split data to a new structure array
34+
hydro_split(i) = temp_hydro;
35+
end
36+
37+
%% interpolate BEM data to greater directional fidelity
38+
theta = wrapTo180(theta);
39+
nTheta0 = length(theta);
40+
thetaInds = 1:nTheta0;
41+
42+
newDirs = -40:0.05:40;
43+
newDirs = setdiff(newDirs,theta); % remove values repeated in theta
44+
45+
% Append the interpolated direction and hydro structue to theta and
46+
% hydro_split respectively.
47+
theta(end+1:end+length(newDirs)) = newDirs;
48+
for i = nTheta0 + 1 : length(theta)
49+
ind1 = thetaInds(theta(i) > theta(1:nTheta0));
50+
ind1 = ind1(end);
51+
52+
ind2 = thetaInds(theta(i) < theta(1:nTheta0));
53+
ind2 = ind2(1);
54+
55+
hydro_split(i) = hydro_split(1);
56+
hydro_split(i).theta = 10; % this has to be the same as the wave direction so that the BEM data processes correctly.
57+
hydro_split(i).file = [hydro.file '_' num2str(wrapTo360(theta(i)))];
58+
59+
dTheta = (theta(i) - theta(ind1)) / (theta(ind2) - theta(ind1));
60+
for iVar = 2:length(vars) % start at 2 to skip theta
61+
hydro_split(i).(vars{iVar}) = hydro_split(ind1).(vars{iVar}) * (1-dTheta) +...
62+
hydro_split(ind2).(vars{iVar}) * dTheta;
63+
end
64+
end
65+
66+
% Sort theta and hydro_split into the correct order based on frequency
67+
[thetaSorted,iSorted] = sort(theta);
68+
theta = wrapTo360(theta);
69+
thetaSorted = wrapTo360(thetaSorted);
70+
hydro_sorted = hydro_split(iSorted);
71+
72+
%% Write all data to h5 files
73+
for i = 1:length(hydro_sorted)
74+
% Skip files that have already been written because writeBEMIOH5 is slow
75+
if ~isfile([hydro_sorted(i).file '.h5'])
76+
writeBEMIOH5(hydro_sorted(i));
77+
end
78+
end

‎Variable_Hydro/Passive_Yaw/irr_passiveYaw.mat‎

2.08 MB
Binary file not shown.

‎Variable_Hydro/Passive_Yaw/plotCases.m‎

Lines changed: 52 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,52 @@
1+
% NOTE: passive yaw and variable hydro forces will not match in surge,
2+
% sway, roll, or pitch because they are being defined in the body's local
3+
% coordinate system, not the global one. This processing could be added
4+
% during the simulation, but is not generalizable to states that are not
5+
% connected to a dof. It's also not necessary to show the current
6+
% comparison to passive yaw and variable hydro.
7+
8+
% To visualize results without runCases.m or re-running simulations, load
9+
% output_regX.mat and reg_passiveYaw.m (or the irregular wave outputs
10+
% respectively)
11+
12+
% Case independent parameters
13+
time = output.bodies(1).time;
14+
elevation = output.wave.elevation;
15+
waves.direction = 10;
16+
17+
% Passive yaw relative position vs BEM direction utilized
18+
figure()
19+
plot(yaw.time, 10-yaw.position*180/pi, yaw.time, yaw.instantDirection, '--');
20+
xlabel('Time (s)');
21+
ylabel('Yaw angle (deg)');
22+
legend('Relative position','Instantaneous BEM direction used');
23+
title('Comparison of relative yaw angle vs BEM data used');
24+
25+
% Compare position and force at various VH discretizations
26+
figure()
27+
t = tiledlayout(1,2);
28+
title(t, 'Passive Yaw vs Variable Hydro Comparison');
29+
xlabel(t, 'Time(s)');
30+
31+
nexttile
32+
plot(yaw.time, yaw.position*180/pi, 'k',...
33+
time, position*180/pi, '--');
34+
grid on
35+
ylabel('Yaw position (deg)')
36+
legend(legendStrings)
37+
38+
nexttile
39+
plot(yaw.time, yaw.forceExcitation, 'k',...
40+
time, forceExcitation, '--');
41+
grid on
42+
ylabel('ExcitationForce in Yaw (N)')
43+
legend(legendStrings)
44+
45+
% Computational time
46+
figure()
47+
bar([yaw.compTime, compTime]/yaw.compTime);
48+
xticklabels(legendStrings);
49+
xlabel('Case');
50+
ylabel('CPU Time / Yaw CPU Time (-), via MATLAB `cputime` function');
51+
title('Computational time comparison');
52+

‎Variable_Hydro/Passive_Yaw/reg_passiveYaw.mat‎

1.74 MB
Binary file not shown.

‎Variable_Hydro/Passive_Yaw/runCases.m‎

Lines changed: 46 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,46 @@
1+
% This script runs a variety of variable hydro resolutions and compares
2+
% outputs to the passive yaw implementation.
3+
4+
% Choose a wave type. This is used by the wecSimInputFile.m
5+
waveFlag = 'regular';
6+
% waveFlag = 'irregular';
7+
8+
if isequal(waveFlag,'regular')
9+
saveStr = 'reg';
10+
elseif isequal(waveFlag,'irregular')
11+
saveStr = 'irr';
12+
end
13+
14+
dThetas = [2, 1, 0.5, 0.25, 0.1, 0.05]; % directional discretization (deg) of the BEM datasets
15+
16+
legendStrings = {'Passive Yaw'};
17+
compTime = zeros(1,length(dThetas));
18+
for itheta = 1:length(dThetas)
19+
% Pre-processing - Iterate through each resolution
20+
dTheta = dThetas(itheta);
21+
legendStrings{itheta + 1} = num2str(dTheta);
22+
23+
bemDirections = -40:dTheta:40;
24+
bemDirections = sort(unique(bemDirections));
25+
26+
% Call WEC-Sim
27+
timeTemp = cputime;
28+
wecSim
29+
compTime(itheta) = cputime - timeTemp;
30+
31+
% Post-processing - save necessary data
32+
hydroForceIndex(itheta,:) = output.bodies(1).hydroForceIndex;
33+
instantDirection(itheta,:) = bemDirections(output.bodies(1).hydroForceIndex);
34+
forceTotal(itheta,:) = output.bodies(1).forceTotal(:,6);
35+
forceExcitation(itheta,:) = output.bodies(1).forceExcitation(:,6);
36+
position(itheta,:) = output.bodies(1).position(:,6);
37+
38+
clear body
39+
save(['output_' saveStr num2str(itheta) '.mat']);
40+
end
41+
42+
% Load in high resolution passive yaw data
43+
load([saveStr '_passiveYaw.mat']);
44+
45+
%% Plots
46+
plotCases

0 commit comments

Comments
 (0)