Add variable mass example

.gitignore

@@ -95,6 +95,8 @@ WECCCOMP/**/SS*.fig
 Passive_Yaw/**/IrrYaw.mat
 Passive_Yaw/**/figYawIrr.fig
 yaw*.mat
+!/Variable_Hydro/Variable_Mass/hydroData/WAMIT/*.h5
+
 
 
 

Variable_Hydro/Variable_Mass/calcIndex.m

@@ -0,0 +1,16 @@
+function body1_hydroForceIndex  = calcIndex(time, hfIndexInitial)
+
+% draftInitial = draftVals(hfIndexInitial);
+% draftInstantaneous = -(zDisp - zCG - draftInitial);
+% [~, body1_hydroForceIndex] = min(abs(draftVals - draftInstantaneous));
+
+% calculate the index based on the current wave conditions
+body1_hydroForceIndex = floor(time/100) + 1; % switch at every 100 seconds
+
+if body1_hydroForceIndex > 9
+    body1_hydroForceIndex = 9;
+end
+
+% body1_hydroForceIndex = 5;
+
+end

Variable_Hydro/Variable_Mass/hydroData/WAMIT/bemio.m

@@ -0,0 +1,17 @@
+close all
+clear all
+clc
+
+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,62 @@
+% This script identifies the dynamics of the float in the respective wave 
+% conditions and determines the optimal proportional gain value for a 
+% passive controller (for regular waves)
+clear all
+close all
+clc
+
+draftVals = 1:9;
+figure()
+
+for ii = 1:length(draftVals)
+
+    % Load hydrodynamic data for float from BEM
+    filename = ['hydroData/WAMIT/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,56 @@
+close all
+
+% Plot heave response and forces for body 1
+output.plotResponse(1,3);
+output.plotForces(1,3);
+
+% find peaks
+% [pks, pksInds] = findpeaks(output.bodies(1).position(:,3));
+% 
+% figure()
+% plot(output.bodies(1).time, output.bodies(1).position(:,3))
+% hold on
+% plot(output.bodies(1).time(pksInds), output.bodies(1).position(pksInds,3),'*')
+% xlabel('time (s)')
+% ylabel('heave disp (m)')
+% 
+% % use 2nd peak 
+% natPeriod = output.bodies(1).time(pksInds(2))/2;
+
+figure()
+plot(hfIndex)
+xlabel('time (s)')
+ylabel('hf index')
+
+figure()
+plot(hfMass)
+xlabel('time (s)')
+ylabel('mass (kg)')
+
+figure()
+plot(output.bodies(1).time, output.bodies(1).position(:,3))
+xlabel('time (s)')
+ylabel('heave disp (m)')
+
+figure()
+plot(output.bodies(1).time,output.bodies(1).forceRestoring(:,3))
+ylabel('restoring force')
+
+figure()
+plot(cgDisp)
+ylabel('disp cg')
+
+figure()
+plot(hsMass)
+ylabel('hs mass')
+
+figure()
+plot(hsVolume)
+ylabel('hs volume')
+
+figure()
+plot(output.ptos(1).time, output.ptos(1).powerInternalMechanics(:,3))
+yline(mean(output.ptos(1).powerInternalMechanics(:,3)))
+text(10,-5e4,sprintf('mean power = %.0f W', mean(output.ptos(1).powerInternalMechanics(:,3))))
+xlabel('time (s)')
+ylabel('power (W)')

Variable_Hydro/Variable_Mass/wecSimInputFile.m

@@ -0,0 +1,64 @@
+%% 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.05;                                 % Simulation Time-Step [s]
+simu.cicEndTime = 15;                           % Specify CI Time [s]
+simu.dtOut = 0.1;
+simu.rho = 1025;
+
+%% Wave Information
+% % noWaveCIC, no waves with radiation CIC  
+% waves = waveClass('noWaveCIC');       % Initialize Wave Class and Specify Type 
+
+% % Regular Waves  
+waves = waveClass('regular');           % Initialize Wave Class and Specify Type                                 
+waves.height = 1;                     % Wave Height [m]
+waves.period = 8;                       % Wave Period [s]
+
+% Waves with imported wave elevation time-history  
+% waves = waveClass('elevationImport');          % Create the Wave Variable and Specify Type
+% waves.elevationFile = 'elevationData.mat';     % Name of User-Defined Time-Series File [:,2] = [time, eta]
+% waves.direction = 0;
+
+%% Body Data
+% Define h5 files for the sphere
+r = 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/WAMIT/draft' num2str(draftVals(ii)), '.h5'];
+
+    fullVolume = (4/3)*pi*r^3;
+    immersedVolume = (1/3)*pi*draftVals(ii)^2*(3*r - draftVals(ii));
+    massVal(ii) = rho*immersedVolume;
+    inertiaVal(ii,:) = (immersedVolume/fullVolume)*2*[20907301 21306090.66 37085481.11];
+end
+
+% sphere
+body(1) = bodyClass(h5Files);  % Initialize bodyClass for Flap
+body(1).geometryFile = 'geometry/sphere.stl';  % Geometry File
+body(1).mass = 'equilibrium';                           % User-Defined mass [kg]
+% body(1).mass = [1e5, 2e5];
+body(1).inertia = inertiaVal(5,:);   % Moment of Inertia [kg*m^2]
+body(1).initial.displacement = [0, 0, 0];
+body(1).variableHydro.option = 1;
+body(1).variableHydro.hydroForceIndexInitial = 1;
+body(1).variableHydro.mass = massVal; % 'equilibrium'
+body(1).variableHydro.inertia = inertiaVal;
+
+% body(1).variableHydro.option = 1;
+% body(1).variableHydro.hydroForceIndexInitial = find(bemDirections==10); % default = 10 deg incident wave
+
+%% 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]