Skip to content

Commit 6df131e

Browse files
TianyuanWangiakeeste
TianyuanWangi
and
akeeste
authored
Radiation damping load extrapolation (WEC-Sim#1476)
* Radiation damping load extrapolation * update convolutionIntegralInterp rename elsewhere * update extrapFrad docstring * update dimensions in docstring * Simpler implementation of the extrapolation feature * Simpler implementation of the extrapolation feature * add extrapFrad to flex body, and library to R2020b * update cicTest format * update size of force output in CIC functions --------- Co-authored-by: akeeste <akeeste@sandia.gov>
1 parent f6595e2 commit 6df131e

6 files changed

Lines changed: 80 additions & 16 deletions

File tree

‎source/functions/initializeWecSim.m‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,7 @@
2424
% Clear old input, plots, log file and start new log file.
2525
diary off
2626
clear body waves simu output pto constraint ptoSim mooring values names InParam
27-
clear ConvolutionIntegral_interp convolutionIntegralSurface % reset functions with persistent variables
27+
clear convolutionIntegralInterp convolutionIntegralSurface % reset functions with persistent variables
2828
try delete('*.log'); end
2929
diary('simulation.log')
3030

source/functions/simulink/model/ConvolutionIntegral_interp.m renamed to source/functions/simulink/model/convolutionIntegralInterp.m

Lines changed: 21 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
function Frad = ConvolutionIntegral_interp(velocity, irkbInput, cicTime)
1+
function [timeExtrap, FradExtrap] = convolutionIntegralInterp(velocity, irkbInput, cicTime, time)
22
%#codegen
33
% Function to calculate convolution integral. velocity is the only dynamic input.
44
% irkb, nDOF and cicTime do not change with time.
@@ -8,7 +8,7 @@
88
% LDOF = radiating dofs from all bodies (6*Nbodies)
99
% nt = length of cicTime (simu.cicEndTime / simu.cicDt)
1010
%
11-
% Paramters:
11+
% Parameters:
1212
% velocity : float [1 LDOF]
1313
% The current velocities of all bodies
1414
% e.g. 6 for 1 body, 12 for 2 bodies and B2B on
@@ -18,21 +18,28 @@
1818
%
1919
% cicTime : float [1 nt]
2020
% All CI times
21+
%
22+
% time : float [1 1]
23+
% The current timestep
2124
%
2225
% Returns:
23-
% Frad : float [nDOF 1]
24-
% Radiation force in each degree of freedom
26+
% timeExtrap : float [3 1]
27+
% Previous 3 main time steps used for extrapolation
28+
% FradExtrap : float [3 nDOF]
29+
% Radiation force in each degree of freedom of the previous 3 main time steps, used for extrapolation
2530
%
2631

2732
% define persistent variables to track velocity_history over time. irkb is
2833
% persistent so that the permuted value is only calculated once.
29-
persistent velocityHistory irkb;
34+
persistent velocityHistory irkb timeHistory FradHistory;
3035

3136
% If this is the first time step (i.e. velocity_history is empty),
3237
% define the persistent variables.
3338
if isempty(velocityHistory)
3439
velocityHistory = zeros(length(cicTime), length(velocity)); % [nt LDOF]
3540
irkb = permute(irkbInput, [1 3 2]); % from [nt nDOF LDOF] to [nt LDOF nDOF]
41+
timeHistory = [-1*(cicTime(2)-cicTime(1));-2*(cicTime(2)-cicTime(1));-3*(cicTime(2)-cicTime(1))];
42+
FradHistory = zeros(3, size(irkbInput, 2));
3643
end
3744

3845
% shift velocity_history and set the first column as the current velocity
@@ -46,6 +53,14 @@
4653
timeSeriesSum = squeeze(sum(timeSeries, 2)); % [nt nDOF]
4754

4855
% integrate over time to get the wave radiation force
49-
Frad = squeeze(trapz(cicTime, timeSeriesSum, 1)); % [nDOF 1]
56+
Frad = squeeze(trapz(cicTime, timeSeriesSum, 1)); % [nDOF 1]
57+
58+
% Prepare the variables used for extrapolation
59+
timeHistory = circshift(timeHistory, 1, 1); % [3 1]
60+
timeHistory(1,:) = time;
61+
FradHistory = circshift(FradHistory, 1, 1); % [3 nDOF]
62+
FradHistory(1,:) = Frad(:)';
63+
timeExtrap = timeHistory;
64+
FradExtrap = FradHistory;
5065

5166
end

‎source/functions/simulink/model/convolutionIntegralSurface.m‎

Lines changed: 20 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
function Frad = convolutionIntegralSurface(velocity, hydroForceIndex, hydroForceIndexInitial, irkbSurfaceInput, cicTime)
1+
function [timeExtrap, FradExtrap] = convolutionIntegralSurface(velocity, hydroForceIndex, hydroForceIndexInitial, irkbSurfaceInput, cicTime, time)
22
%#codegen
33
% Function to calculate convolution integral from a surface varying in
44
% time, DOF, and variable hydro state.
@@ -26,16 +26,21 @@
2626
%
2727
% cicTime : float [1 nt]
2828
% All CI times
29+
%
30+
% time : float [1 1]
31+
% The current timestep
2932
%
3033
% Returns:
31-
% Frad : float [nDOF 1]
32-
% Radiation force in each degree of freedom
34+
% timeExtrap : float [3 1]
35+
% Previous 3 main time steps used for extrapolation
36+
% FradExtrap : float [3 nDOF]
37+
% Radiation force in each degree of freedom of the previous 3 main time steps, used for extrapolation
3338
%
3439

3540
% define persistent variables to track velocity_history and
3641
% hydroForceIndexSurface over time. irkb is persistent so that the permuted
3742
% value is only calculated once.
38-
persistent velocityHistory irkbSurface hydroForceIndexSurface;
43+
persistent velocityHistory irkbSurface hydroForceIndexSurface timeHistory FradHistory;
3944

4045
% If this is the first time step (i.e. velocity_history is empty),
4146
% define the persistent variables.
@@ -47,6 +52,8 @@
4752
for i = 1:size(hydroForceIndexSurface, 1)
4853
hydroForceIndexSurface(i, 1, 1, hydroForceIndexInitial) = true;
4954
end
55+
timeHistory = [-1*(cicTime(2)-cicTime(1));-2*(cicTime(2)-cicTime(1));-3*(cicTime(2)-cicTime(1))];
56+
FradHistory = zeros(3, size(irkbSurfaceInput, 2));
5057
end
5158

5259
% shift velocity_history and set the first column as the current velocity
@@ -68,6 +75,14 @@
6875
timeSeriesSum = squeeze(sum(timeSeries, 2)); % [nt nDOF]
6976

7077
% integrate over time to get the wave radiation force
71-
Frad = squeeze(trapz(cicTime, timeSeriesSum, 1)); % [nDOF 1]
78+
Frad = squeeze(trapz(cicTime, timeSeriesSum, 1)); % [nDOF 1]
79+
80+
% Prepare the variables used for extrapolation
81+
timeHistory = circshift(timeHistory, 1, 1); % [3 1]
82+
timeHistory(1,:) = time;
83+
FradHistory = circshift(FradHistory, 1, 1); % [3 nDOF]
84+
FradHistory(1,:) = Frad(:)';
85+
timeExtrap = timeHistory;
86+
FradExtrap = FradHistory;
7287

7388
end

‎source/functions/simulink/model/extrapFrad.m‎

Lines changed: 31 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,31 @@
1+
function Frad = extrapFrad(timeExtrap, FradExtrap, time)
2+
%#codegen
3+
% This function is to extrapolate the radiation damping loads when using CIC.
4+
% Since CIC uses persistent variables to track velocityHistory over time,
5+
% radiation damping loads cannot be updated at each sub-time step
6+
%
7+
% This function must be separated from Simulink calls to
8+
% `convolutionIntegralInterp` and `convolutionIntegralSurface` so that its
9+
% block's timestep can be inherited (-1) and be called at minor time steps.
10+
%
11+
% Dimensions:
12+
% nDOF = the body's number of degrees of freedom = body.dof
13+
%
14+
% Parameters:
15+
% timeExtrap : float [3 1]
16+
% Previous 3 main time steps used for extrapolation
17+
%
18+
% FradExtrap : float [3 nDOF]
19+
% Radiation force in each degree of freedom of the previous 3 main time steps, used for extrapolation
20+
%
21+
% time : float [1 1]
22+
% The current time step
23+
%
24+
% Returns:
25+
% Frad : float [nDOF 1]
26+
% The radiation force in each degree of freedom, extrapolated to the current time step
27+
%
28+
29+
Frad = interp1(timeExtrap,FradExtrap,time,'pchip','extrap')';
30+
31+
end

‎source/lib/WEC-Sim/WECSim_Lib_Body_Elements.slx‎

16 KB
Binary file not shown.

‎tests/cicTest.m‎

Lines changed: 7 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -51,10 +51,11 @@
5151

5252
methods(TestClassSetup)
5353
function calcForceRad(testCase)
54-
clear ConvolutionIntegral_interp
54+
clear convolutionIntegralInterp
5555
tmp = cputime;
5656
for it = 1:size(testCase.velocity,1)
57-
testCase.Frad0(it,:) = ConvolutionIntegral_interp(testCase.velocity(it,:), testCase.irkbSurfaceInput(:,:,:,1), testCase.cicTime);
57+
[~,fRad] = convolutionIntegralInterp(testCase.velocity(it,:), testCase.irkbSurfaceInput(:,:,:,1), testCase.cicTime, testCase.time(it));
58+
testCase.Frad0(it,:) = fRad(1);
5859
end
5960
testCase.t0 = cputime - tmp;
6061
end
@@ -64,15 +65,17 @@ function calcForceRadSurface(testCase)
6465
clear convolutionIntegralSurface
6566
tmp = cputime;
6667
for it = 1:size(testCase.velocity,1)
67-
testCase.FradS1(it,:) = convolutionIntegralSurface(testCase.velocity(it,:), 1, 1, testCase.irkbSurfaceInput, testCase.cicTime);
68+
[~,fRad] = convolutionIntegralSurface(testCase.velocity(it,:), 1, 1, testCase.irkbSurfaceInput, testCase.cicTime, testCase.time(it));
69+
testCase.FradS1(it,:) = fRad(1);
6870
end
6971
testCase.ts1 = cputime - tmp;
7072

7173
% Test CI surface calculation while switching variable hydro states
7274
clear convolutionIntegralSurface
7375
tmp = cputime;
7476
for it = 1:size(testCase.velocity,1)
75-
testCase.FradS2(it,:) = convolutionIntegralSurface(testCase.velocity(it,:), testCase.hydroForceIndex(it), 1, testCase.irkbSurfaceInput, testCase.cicTime);
77+
[~,fRad] = convolutionIntegralSurface(testCase.velocity(it,:), testCase.hydroForceIndex(it), 1, testCase.irkbSurfaceInput, testCase.cicTime, testCase.time(it));
78+
testCase.FradS2(it,:) = fRad(1);
7679
end
7780
testCase.ts2 = cputime - tmp;
7881
end

0 commit comments

Comments
 (0)