User defined force - RM3

RM3_UserDefinedForce/PlotResults.m

@@ -0,0 +1,17 @@
+data3 = load("ResultsV1\output_force3.mat");
+data0 = load("output_noForce.mat");
+
+figure()
+plot(data0.output.bodies(1).time, data0.output.bodies(1).position(:,4),output.bodies(1).time, output.bodies(1).position(:,4),data3.output.bodies(1).time, data3.output.bodies(1).position(:,4))
+legend('No Force', 'Ramp', 'No Ramp')
+title('Roll - X rotation')
+
+figure()
+plot(data0.output.bodies(1).time, data0.output.bodies(1).position(:,5),output.bodies(1).time, output.bodies(1).position(:,5),data3.output.bodies(1).time, data3.output.bodies(1).position(:,5))
+legend('No Force', 'Ramp', 'No Ramp')
+title('Pitch - Y rotation')
+
+figure()
+plot(data0.output.bodies(1).time, data0.output.bodies(1).position(:,6),output.bodies(1).time, output.bodies(1).position(:,6),data3.output.bodies(1).time, data3.output.bodies(1).position(:,6))
+legend('No Force', 'Ramp', 'No Ramp')
+title('Yaw - Z rotation')

RM3_UserDefinedForce/README.md

@@ -0,0 +1,18 @@
+# RM3 User-Defined Force
+
+**Author:** Jorge Leon & Rebekah A. Saucier
+
+**Geometry:** RM3
+
+**Original Version:** WEC-Sim v7.1
+
+**Description:**
+
+The RM3 example was modified to include an external force defined by the user.
+
+This example demonstrates how a user-defined external force can be applied body 1 of to the RM3 model in WEC-Sim. The external force is defined in `wecSimInputFile.m`  as:
+
+`F_ext_b1`
+
+
+

RM3_UserDefinedForce/TestRM3UserDefinedForce.m

@@ -0,0 +1,149 @@
+classdef TestRM3UserDefinedForce < matlab.unittest.TestCase
+
+    properties
+        OriginalDefault
+        testDir
+    end
+
+    methods (Access = 'public')
+        function obj = TestRM3UserDefinedForce
+            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
+    end
+
+    methods (TestClassTeardown)
+        function checkVisibilityRestored(testCase)
+            set(0, 'DefaultFigureVisible', testCase.OriginalDefault);
+            testCase.assertEqual(get(0, 'DefaultFigureVisible'), ...
+                                 testCase.OriginalDefault);
+        end
+    end
+
+    methods (Test)
+        function testRM3UserDefinedForceRuns(testCase)
+            currentDir = pwd;
+            cleanup = onCleanup(@() cd(currentDir));
+
+            cd(testCase.testDir);
+
+            wecSim
+        end
+
+        function testUserDefinedForceInput(testCase)
+            F_ext_b1 = testCase.loadUserDefinedForce();
+
+            % F_ext_b1 should be a generalized force/moment vector:
+            % [Fx, Fy, Fz, Mx, My, Mz]
+            testCase.verifyTrue(isnumeric(F_ext_b1), ...
+                'F_ext_b1 should be numeric.');
+
+            testCase.verifySize(F_ext_b1, [1 6], ...
+                'F_ext_b1 should be a 1-by-6 force/moment vector.');
+
+            testCase.verifyTrue(all(isfinite(F_ext_b1(:))), ...
+                'F_ext_b1 should contain only finite values.');
+
+            testCase.verifyTrue(norm(F_ext_b1) > 0, ...
+                'F_ext_b1 should define a nonzero external force or moment.');
+        end
+
+        function testEquivalentForceMomentCouple(testCase)
+            F_ext_b1 = testCase.loadUserDefinedForce();
+
+            % Split generalized force into translational force and moment.
+            % F_ext_b1 = [Fx, Fy, Fz, Mx, My, Mz]
+            F = F_ext_b1(1:3).';
+            M = F_ext_b1(4:6).';
+
+            testCase.verifyTrue(norm(F) > 0, ...
+                'F_ext_b1 should include a nonzero translational force.');
+
+            applicationPoints = testCase.getApplicationPoints();
+
+            for i = 1:size(applicationPoints, 1)
+                % r_CG_P is the vector from the body CG to the force
+                % application point P.
+                r_CG_P = applicationPoints(i, :).';
+
+                % Force applied at P, represented as an equivalent wrench
+                % about the body CG.
+                wrenchFromForceAtPoint = ...
+                    testCase.forceAtPointToWrenchAtCG(F, M, r_CG_P);
+
+                % Same force applied at the body CG, plus equivalent moment
+                % r_CG_P x F.
+                equivalentForceMomentCoupleAtCG = [
+                    F
+                    M + cross(r_CG_P, F)
+                ];
+
+                testCase.verifyEqual(wrenchFromForceAtPoint, ...
+                                     equivalentForceMomentCoupleAtCG, ...
+                                     'AbsTol', 1e-12);
+            end
+        end
+    end
+
+    methods (Access = private)
+        function F_ext_b1 = loadUserDefinedForce(testCase)
+            currentDir = pwd;
+            cleanup = onCleanup(@() cd(currentDir));
+
+            cd(testCase.testDir);
+
+            run('wecSimInputFile.m');
+
+            testCase.assertTrue(exist('F_ext_b1', 'var') == 1, ...
+                'F_ext_b1 should be defined in wecSimInputFile.m.');
+        end
+
+        function applicationPoints = getApplicationPoints(~)
+            % Representative application points in the model Cartesian
+            % coordinate system. Each row is [x, y, z] in meters and
+            % represents an offset from the body CG to force application
+            % point P.
+            %
+            % These are not RM3-specific geometry points. They are used to
+            % verify the equivalent force-moment relationship for multiple
+            % possible offsets.
+            applicationPoints = [
+                 0.0,  0.0,  0.0
+                 1.0,  0.0,  0.0
+                 0.0,  1.0,  0.0
+                 0.0,  0.0,  1.0
+                -1.0,  0.0,  0.0
+                 0.0, -1.0,  0.0
+                 0.0,  0.0, -1.0
+                 1.0,  1.0,  0.0
+                 1.0,  0.0,  1.0
+                 0.0,  1.0,  1.0
+                 2.0, -1.0,  0.5
+                -0.5,  3.0, -1.0
+            ];
+        end
+
+        function wrench = forceAtPointToWrenchAtCG(~, F, M, r_CG_P)
+            % Equivalent wrench at the body CG for a force F applied at
+            % point P:
+            %
+            % M_CG = M + r_CG_P x F
+
+            wrench = [
+                F
+                M + cross(r_CG_P, F)
+            ];
+        end
+    end
+end

RM3_UserDefinedForce/dataVisualization.m

@@ -0,0 +1,25 @@
+%% Data Visualization: Force Comparison (Updated Structure)
+% Load the datasets
+data0 = load("output_noForce.mat");
+data1 = load("output_force1.mat");
+data2 = load("output_force2.mat");
+data3 = load("output_force3.mat");
+data4 = load("output_force4.mat");
+
+figure('Color', 'w'); 
+hold on; grid on;
+
+% Plot each line using the updated .output.bodies path
+plot(data0.output.bodies(1).time, data0.output.bodies(1).position(:,3), 'k--', 'LineWidth', 1.5);
+plot(data1.output.bodies(1).time, data1.output.bodies(1).position(:,3), 'LineWidth', 1.2);
+plot(data2.output.bodies(1).time, data2.output.bodies(1).position(:,3), 'LineWidth', 1.2);
+plot(data3.output.bodies(1).time, data3.output.bodies(1).position(:,3), 'LineWidth', 1.2);
+plot(data4.output.bodies(1).time, data4.output.bodies(1).position(:,3), 'LineWidth', 1.2);
+
+% Formatting the plot
+xlabel('Time (s)', 'FontSize', 12);
+ylabel('Vertical Position (m)', 'FontSize', 12);
+title('Body 1 Heave Response Across Force Profiles', 'FontSize', 14);
+
+legend('No Force', 'Force 1', 'Force 2', 'Force 3', 'Force 4', ...
+       'Location', 'best');

RM3_UserDefinedForce/products.txt

@@ -0,0 +1,3 @@
+Simulink
+Simscape
+Simscape_Multibody

RM3_UserDefinedForce/userDefinedFunctions.m

@@ -0,0 +1,25 @@
+%Example of user input MATLAB file for post processing
+
+%Plot waves
+waves.plotElevation(simu.rampTime);
+try 
+    waves.plotSpectrum();
+catch
+end
+
+%Plot heave response for body 1
+output.plotResponse(1,3);
+
+%Plot heave response for body 2
+output.plotResponse(2,3);
+
+%Plot heave forces for body 1
+output.plotForces(1,3);
+
+%Plot heave forces for body 2
+output.plotForces(2,3);
+
+%Save waves and response as video
+% output.saveViz(simu,body,waves,...
+%     'timesPerFrame',5,'axisLimits',[-150 150 -150 150 -50 20],...
+%     'startEndTime',[100 125]);

RM3_UserDefinedForce/wecSimInputFile.m

@@ -0,0 +1,50 @@
+%% Simulation Data
+simu = simulationClass();               % Initialize Simulation Class
+simu.simMechanicsFile = 'RM3.slx';      % Specify Simulink Model File
+simu.mode = 'normal';                   % Specify Simulation Mode ('normal','accelerator','rapid-accelerator')
+simu.explorer = 'off';                   % Turn SimMechanics Explorer (on/off)
+simu.startTime = 0;                     % Simulation Start Time [s]
+simu.rampTime = 100;                    % Wave Ramp Time [s]
+simu.endTime = 400;                     % Simulation End Time [s]
+simu.solver = 'ode4';                   % simu.solver = 'ode4' for fixed step & simu.solver = 'ode45' for variable step 
+simu.dt = 0.1; 							% Simulation time-step [s]
+
+%% 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 = 2.5;                     % Wave Height [m]
+waves.period = 8;                       % Wave Period [s]
+
+%% Body Data
+% Float
+%External force on Body 1
+F_ext_b1 = [0,0,5.0e5,0,0,0];
+
+body(1) = bodyClass('../_Common_Input_Files/RM3/hydroData/rm3.h5');
+    % Create the body(1) Variable, Set Location of Hydrodynamic Data File
+    % and Body Number Within this File.
+body(1).geometryFile = '../_Common_Input_Files/RM3/geometry/float.stl';    % Location of Geometry File
+body(1).mass = 'equilibrium';
+    % Body Mass. The 'equilibrium' Option Sets it to the Displaced Water
+    % Weight.
+body(1).inertia = [20907301 21306090.66 37085481.11];  % Moment of Inertia [kg*m^2]     
+
+% Spar/Plate
+body(2) = bodyClass('../_Common_Input_Files/RM3/hydroData/rm3.h5');
+body(2).geometryFile = '../_Common_Input_Files/RM3/geometry/plate.stl';
+body(2).mass = 'equilibrium';                   
+body(2).inertia = [94419614.57 94407091.24 28542224.82];
+
+%% PTO and Constraint Parameters
+% Floating (3DOF) Joint
+constraint(1) = constraintClass('Constraint1'); % Initialize Constraint Class for Constraint1
+constraint(1).location = [0 0 0];               % Constraint Location [m]
+
+% Translational PTO
+pto(1) = ptoClass('PTO1');                      % Initialize PTO Class for PTO1
+pto(1).stiffness = 0;                           % PTO Stiffness [N/m]
+pto(1).damping = 1200000;                       % PTO Damping [N/(m/s)]
+pto(1).location = [0 0 0];                      % PTO Location [m]