Add variable mass example

Variable_Hydro/Passive_Yaw/TestVariableHydro.m → Variable_Hydro/TestVariableHydro.m

@@ -3,12 +3,14 @@
     properties
         OriginalDefault
         testDir
-        h5Dir = fullfile("hydroData")
-        h5Name = 'oswec_0.h5'
+        h5DirPY = ['Passive_Yaw',filesep,'hydroData']
+        h5NamePY = 'oswec_0.h5'
         outName = 'oswec.out'
+        h5DirVM = ['Variable_Mass',filesep,'hydroData']
+        h5NameVM = 'draft1.h5'
+
     end
 
-
     methods (Access = 'public')
         function obj = TestVariableHydro
             obj.testDir = fileparts(mfilename('fullpath'));
@@ -27,9 +29,19 @@ function captureVisibility(testCase)
             testCase.OriginalDefault = get(0, 'DefaultFigureVisible');
         end
 
-        function runBemio(testCase)
-            cd(testCase.h5Dir);
-            if isfile(testCase.h5Name)
+        function runBemioPY(testCase)
+            cd(testCase.h5DirPY);
+            if isfile(testCase.h5NamePY)
+                fprintf('runBemio skipped, *.h5 already exists\n')
+            else
+                bemio
+            end
+            cd(testCase.testDir)
+        end
+
+        function runBemioVM(testCase)
+            cd(testCase.h5DirVM);
+            if isfile(testCase.h5NameVM)
                 fprintf('runBemio skipped, *.h5 already exists\n')
             else
                 bemio
@@ -38,21 +50,30 @@ function runBemio(testCase)
         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 
     end
 
     methods(Test)
-        function testCable(testCase)
+        function testPY(testCase)
+            cd('Passive_Yaw')
             runCases
         end
+        function testVM(testCase)
+            cd('Variable_Mass')
+            wecSim
+        end
     end
 
 end

Variable_Hydro/Variable_Mass/README.md

@@ -0,0 +1,20 @@
+# Variable Hydrodynamics Variable Mass
+
+**Author:** Jeff Grasberger
+
+**Geometry:**	Sphere
+
+**Original Version:** 	WEC-Sim v6.1
+
+**Description**
+
+This application demonstrates using WEC-Sim's variable hydrodynamics feature to model a sphere with a variable mass and draft.
+The sphere's mass is set-up to change instantaneously every 100 seconds, increasing the its draft. 
+The existing model is simple but can be built upon to demonstrate any WEC that changes mass during a simulation.
+WEC-Sim by default does not support variable mass, so this model requires a broken library link to add the General Variable Mass block.
+To recreate (for future library updates): duplicate the changes within sphereVarMass/Float/Hydrodynamic Body/Structure
+R2021b or later is required because the Simulink model uses the "Create Diagonal Matrix" block which was only moved to standard Simulink in R2021b.
+
+**Relevant Citation(s)**
+
+Keester, A.; Ogden, D.; Husain, S.; Ruehl, K.; Pan, J.; Quiroga, J.; Grasberger, J.; Forbush, D.; Tom, N.; Housner, S.; Tran, T. (2023). Review of TEAMER Awards for WEC-Sim Support. Paper presented at 15th European Wave and Tidal Energy Conference (EWTEC 2023), Bilbao, Spain. https://doi.org/10.36688/ewtec-2023-261

Variable_Hydro/Variable_Mass/calcIndex.m

@@ -0,0 +1,9 @@
+function body1_hydroForceIndex = calcIndex(time)
+
+body1_hydroForceIndex = floor(time/100) + 1; % Increase the mass every 100 seconds
+
+if body1_hydroForceIndex > 9
+    body1_hydroForceIndex = 9; % restrict to 9 (# of hydrodata files)
+end
+
+end

Variable_Hydro/Variable_Mass/hydroData/bemio.m

@@ -0,0 +1,15 @@
+close all
+
+for ii = 1:9
+
+    hydro = struct();
+
+    outName = ['draft' num2str(ii) '.out'];
+
+    hydro = readWAMIT(hydro,outName,[]);
+    hydro = radiationIRF(hydro,15,[],[],[],6);
+    hydro = radiationIRFSS(hydro,[],[]);
+    hydro = excitationIRF(hydro,15,[],[],[],[]);
+    writeBEMIOH5(hydro)
+
+end

Variable_Hydro/Variable_Mass/impedanceAnalysis.m

@@ -0,0 +1,57 @@
+% This script identifies the dynamics of the float for the different draft
+% values. This script is not used by the example here but could be used to 
+% design a controller that changes the draft for different wave conditions.
+draftVals = 1:9;
+figure()
+
+for ii = 1:length(draftVals)
+    % Load hydrodynamic data for float from BEM
+    filename = ['hydroData/draft' num2str(draftVals(ii)) '.h5'];
+    hydro = readH5ToStruct(filename);
+    rho = 1000;
+    gravity = 9.81;
+
+    ptoDamping = 0;
+
+    % Define the intrinsic mechanical impedance for the device
+    mass = rho*hydro.Vo;
+    addedMass = squeeze(hydro.A(3,3,:))*rho;
+    radiationDamping = squeeze(hydro.B(3,3,:)).*squeeze(hydro.w')*rho + ptoDamping;
+    hydrostaticStiffness = hydro.Khs(3,3)*rho*gravity;
+    Gi = -((hydro.w)'.^2.*(mass+addedMass)) + 1j*hydro.w'.*radiationDamping + hydrostaticStiffness;
+    Zi = Gi./(1j*hydro.w');
+
+    % Calculate magnitude and phase for bode plot
+    Mag = 20*log10(abs(Zi));
+    Phase = (angle(Zi))*(180/pi);
+
+    % calculate natural frequency
+    [~,closestIndNat] = min(abs(imag(Zi)));
+    natFreqs(ii) = hydro.w(closestIndNat)/(2*pi); % Hz
+
+    % Create bode plot for impedance
+    subplot(2,1,1)
+    semilogx((hydro.w)/(2*pi),Mag)
+    hold on
+
+    subplot(2,1,2)
+    semilogx((hydro.w)/(2*pi),Phase)
+    hold on
+
+    legendLabel{ii} = ['draft = ' num2str(draftVals(ii)) ' m'];
+end
+
+subplot(2,1,1)
+xlabel('freq (hz)')
+ylabel('mag (dB)')
+xlim([1e-2, 1e0])
+ylim([80, 140])
+grid on
+
+subplot(2,1,2)
+xlabel('freq (hz)')
+ylabel('phase (deg)')
+grid on
+xlim([1e-2, 1e0])
+ylim([-100, 100])
+legend(legendLabel,"Location","west")

Variable_Hydro/Variable_Mass/userDefinedFunctions.m

@@ -0,0 +1,27 @@
+close all
+
+% Plot heave response and forces for body 1
+output.plotResponse(1,3);
+output.plotForces(1,3);
+
+figure()
+plot(output.bodies(1).time, output.bodies(1).hydroForceIndex)
+xlabel('Time (s)')
+ylabel('Index (-)');
+title('Sphere hydroForceIndex');
+
+hfNames = fieldnames(body.hydroForce);
+for ii = 1:length(hfNames)
+    massVec(ii) = body.hydroForce.(hfNames{ii}).mass;
+end
+
+figure()
+plot(output.bodies(1).time, massVec(output.bodies(1).hydroForceIndex))
+xlabel('Time (s)')
+ylabel('Mass (kg)')
+title('Sphere mass');
+
+figure()
+plot(output.ptos(1).time, output.ptos(1).powerInternalMechanics(:,3))
+xlabel('Time (s)')
+ylabel('Power (W)')

Variable_Hydro/Variable_Mass/wecSimInputFile.m

@@ -0,0 +1,52 @@
+%% Simulation Data
+simu = simulationClass();                       % Initialize Simulation Class
+simu.simMechanicsFile = 'sphereVarMass.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.rampTime = 0;                              % Wave Ramp Time [s]
+simu.endTime = 900;                             % Simulation End Time [s]        
+simu.solver = 'ode4';                           % simu.solver = 'ode4' for fixed step & simu.solver = 'ode45' for variable step 
+simu.dt = 0.01;                                 % Simulation Time-Step [s]
+simu.cicEndTime = 15;                           % Specify CI Time [s]
+simu.dtOut = 0.1;
+simu.rho = 1025;
+
+%% Wave Information
+% Regular Waves  
+waves = waveClass('regular');           % Initialize Wave Class and Specify Type                                 
+waves.height = 1;                     % Wave Height [m]
+waves.period = 8;                       % Wave Period [s]
+
+%% Body Data
+% Define h5 files for the sphere
+radius = 5;
+rho = 1025;
+draftVals = 1:9;
+natFreqs = [0.3979    0.3247    0.2865    0.2546    0.2292    0.2037    0.1751    0.1401    0.0987]; % from impedanceAnalysis.m
+
+for ii = 1:length(draftVals)
+    h5Files{ii} = ['hydroData/draft' num2str(draftVals(ii)), '.h5'];
+
+    fullVolume = (4/3)*pi*radius^3;
+    immersedVolume = (1/3)*pi*draftVals(ii)^2*(3*radius - draftVals(ii));
+    massVals(ii) = rho*immersedVolume; 
+    inertiaVals(ii,:) = (immersedVolume/fullVolume)*2*[20907301 21306090.66 37085481.11];
+end
+
+% Sphere
+body(1) = bodyClass(h5Files);
+body(1).geometryFile = '../../_Common_Input_Files/Sphere/geometry/sphere.stl';
+body(1).mass = 'equilibrium';
+body(1).inertia = inertiaVals(5,:);
+body(1).initial.displacement = [0, 0, 0];
+body(1).variableHydro.option = 1;
+body(1).variableHydro.hydroForceIndexInitial = 1;
+body(1).variableHydro.mass = massVals; % 'equilibrium' for each sphere size. If not set, equilibrium will be assumed by WEC-Sim.
+body(1).variableHydro.inertia = inertiaVals;
+
+%% PTO and Constraint Parameters
+pto(1) = ptoClass('PTO1');                   % Initialize ptoClass for PTO1
+pto(1).stiffness = 0;                        % PTO Stiffness Coeff [Nm/rad]
+pto(1).damping = 2e5;                        % PTO Damping Coeff [Nsm/rad]
+pto(1).location = [0 0 0];                   % PTO Location [m]