diff --git a/RM3_UserDefinedForce/PlotResults.m b/RM3_UserDefinedForce/PlotResults.m new file mode 100644 index 00000000..485230c5 --- /dev/null +++ b/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') \ No newline at end of file diff --git a/RM3_UserDefinedForce/README.md b/RM3_UserDefinedForce/README.md new file mode 100644 index 00000000..12d7999d --- /dev/null +++ b/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` + + + diff --git a/RM3_UserDefinedForce/RM3.slx b/RM3_UserDefinedForce/RM3.slx new file mode 100644 index 00000000..61edfabb Binary files /dev/null and b/RM3_UserDefinedForce/RM3.slx differ diff --git a/RM3_UserDefinedForce/ResultsV1/output_force1.mat b/RM3_UserDefinedForce/ResultsV1/output_force1.mat new file mode 100644 index 00000000..d6f10f54 Binary files /dev/null and b/RM3_UserDefinedForce/ResultsV1/output_force1.mat differ diff --git a/RM3_UserDefinedForce/ResultsV1/output_force2.mat b/RM3_UserDefinedForce/ResultsV1/output_force2.mat new file mode 100644 index 00000000..de5a285a Binary files /dev/null and b/RM3_UserDefinedForce/ResultsV1/output_force2.mat differ diff --git a/RM3_UserDefinedForce/ResultsV1/output_force3.mat b/RM3_UserDefinedForce/ResultsV1/output_force3.mat new file mode 100644 index 00000000..4519263f Binary files /dev/null and b/RM3_UserDefinedForce/ResultsV1/output_force3.mat differ diff --git a/RM3_UserDefinedForce/ResultsV1/output_force4.mat b/RM3_UserDefinedForce/ResultsV1/output_force4.mat new file mode 100644 index 00000000..2377f33a Binary files /dev/null and b/RM3_UserDefinedForce/ResultsV1/output_force4.mat differ diff --git a/RM3_UserDefinedForce/TestRM3UserDefinedForce.m b/RM3_UserDefinedForce/TestRM3UserDefinedForce.m new file mode 100644 index 00000000..3421c622 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/RM3_UserDefinedForce/dataVisualization.m b/RM3_UserDefinedForce/dataVisualization.m new file mode 100644 index 00000000..cbf06b0c --- /dev/null +++ b/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'); \ No newline at end of file diff --git a/RM3_UserDefinedForce/output_noForce.mat b/RM3_UserDefinedForce/output_noForce.mat new file mode 100644 index 00000000..9d034cb8 Binary files /dev/null and b/RM3_UserDefinedForce/output_noForce.mat differ diff --git a/RM3_UserDefinedForce/products.txt b/RM3_UserDefinedForce/products.txt new file mode 100644 index 00000000..82ab62e2 --- /dev/null +++ b/RM3_UserDefinedForce/products.txt @@ -0,0 +1,3 @@ +Simulink +Simscape +Simscape_Multibody diff --git a/RM3_UserDefinedForce/userDefinedFunctions.m b/RM3_UserDefinedForce/userDefinedFunctions.m new file mode 100644 index 00000000..57cc8844 --- /dev/null +++ b/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]); diff --git a/RM3_UserDefinedForce/wecSimInputFile.m b/RM3_UserDefinedForce/wecSimInputFile.m new file mode 100644 index 00000000..893b75ba --- /dev/null +++ b/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]