Asked by Ryan O'Connor
on 23 Apr 2019

I am trying to model a buoy in sinesoidal waves. I am using function handles to do this. There is 4 if conditions:

- When there is no itersection between the buoy and the wave and the buoy is obove water
- When there is no itersection between the buoy and the wave and the buoy is below water
- When there is itersection between the buoy and the wave less that 50%
- When there is no itersection between the buoy and the wave greater than 50%

However what occurs is that it runs through the inital values and determines which if statement to use in the for loop and then just uses that value, even though the values effecting the if statement above the for loop are changing after each loop.

tf_f = @(T,Y) isempty(p_f(T,Y)); % Is the matrix empty == true(1) or false (0) [ is there intersection points from the buoy and wave surface]

%% Conditions for forces

% float

if tf_f(T,Y) == 1 && bottom_f(Y) > p21top_f(T) %if intersection matrix is empty and bottom of buoy is above the wave surface

elseif tf_f(T,Y) == 1 && top_f(Y) < p21top_f(T) %if intersection matrix is empty and top of buoy is below the wave surface

elseif tf_f(T,Y) == 0 && centre_check_f(Y) > p21top_f(T) %if intersection matrix is full and centre of buoy is above the wave surface

elseif tf_f(T,Y) == 0 && centre_check_f(Y) < p21top_f(T) %if intersection matrix is full and centre of buoy is below the wave surface

end

Above is the if statement that decides the how it should act, tf_f is used to determine whether the intersection matrix is empty, however it only does the above for the inital and does not update? after this the Runga_kutta 4th order intergration takesplace in a for loop

%Float (centre)

Accel_h_f = @(T,Y,V) ((Sum_forces_h_f(T,Y,V))/(abs(mass_f)));

for i=1:(length(t)-1) % calculation loop

k_1vh_f = Accel_h_f(t(i), y_f(i), v_f(i))*h;

From what i can tell, all other parts of my code are correct and for each of the conditions above. I just cant seem to figure out how to get the code to check the if statement each time. Any suggestions?

Answer by dpb
on 23 Apr 2019

Edited by dpb
on 24 Apr 2019

tf_f = @(T,Y) isempty(p_f(T,Y));

if tf_f(T,Y) == 1 && bottom_f(Y) > p21top_f(T)

elseif tf_f(T,Y) == 1 && top_f(Y) < p21top_f(T)

elseif tf_f(T,Y) == 0 && centre_check_f(Y) > p21top_f(T)

elseif tf_f(T,Y) == 0 && centre_check_f(Y) < p21top_f(T)

end

Not enough code to be able to tell...we don't know where the loop is which is crucial...

But, the crystal ball is back from the shop so let's give it a whirl and see if it is (finally!) repaired or not--

It says, since the content of any variable not passed as a dummy argument into an anonymous function is the value for that variable in the workspace at the time the anonymous function is defined and is invariant, that definition is outside the loop and so the array p_f is constant in the function definition instead of being revised as you expect.

Either pass p_f, too, or don't use the anonymous function but just write the test expression.

(Now to wait and see if I wasted money on the repair bill... :) )

Ryan O'Connor
on 28 Apr 2019

Apologies for the late repsonse; I was busy with exams. I tried moving the for loop above the if statement and it worked as far as it now changes between the conditions.

However, there is now a seperate issue. While the buoy and the wave are intersecting tf_f = 0 the interx function doesnt appear to find these points causing it to error out. It seems that it fails on the last step in before transitioning, when it is about to transittion from tf=0 to tf = 1.

Index in position 2 exceeds array bounds.

Error in Floating_perfect_Save_Spare>@(Q,Z,S)Q(Z,S)

Error in Floating_perfect_Save_Spare>@(T,Y)select(p(T,Y),1,2)

Although if you look at the plot you can clearly see that there is still intersection so I can't seem to determine why its missing finding these point??

I've tried including a try and catch function. Which to my knowlegde is if it errors do this but it never did the action.

Any idea what could be causing the issue or anyways around it? I've included a good portion of the code so you can see it for yourself since I left out a lot last time

%% Preparation

clc

clear all

close all

%% inital Conditions to be used in the function handles

h=0.005; % step size

t = 0:h:3; % Calculates upto y(3)

y = zeros(1,length(t));

v = zeros(1,length(t));

y(1) = 0;

v(1) = 0;

Y = y(1);

T = t(1);

v = v(1);

%% Constants

g = 9.81; % gravity [m/s2]

DSw = 1025; % Density sea water [kg/m3]

L = 1; % wavelength [m]

A = 0.0375; % wave amplitude [m]

Ti = 1; % Wave Period [s]

c = L/Ti; % wave speed (celerity)[m/s]

d = 0.45; % Water Depth [m]

H = 2*A; %wave height [m]

k = (2*pi)/L; % Wave Numebr

w = (2*pi)/Ti; % Wave Frequency [rads/s]

wl = 0; % water line for wave velocity and acceleration

S_k = 90; % Spring constant [N/m]

%% Coefficents

% Sphere

Cd = 0.1; % Drag

eh = 0.5; %radiation damping heave

es = 0.63; %radiation damping surge

amh = 0.55; % added mass heave

ams = 0.63; % added mass surge

%% Buoy parameters

r=0.075; % radius of buoy

Vbuoy_Total = 4/3*pi*r^3; % total volume of the buoy assuming spherical [m3]

Denpercent = 0.0635;

Den_buoy = DSw*Denpercent;

mass = Den_buoy*Vbuoy_Total; % mass of buoy[kg]

Weight = mass*g; % weight of buoy [N]

Fbmax = DSw*g*Vbuoy_Total; % max possible bupoyancy fully submerged [N]

numwave = 2; % number of desired waves on plot

inc = 150; % Number of increment for step size for time and wavelength

inct = Ti/inc;

ts = 0:inct:numwave*Ti; % time [s]

inclambda = L/inc; % increment step size for a single wave length

ts = 0:inclambda:numwave*L;

th = 0:pi/((numwave*inc)/2):2*pi; % phase angle

%% for loop begin

for i=1:(length(t)-1) % calculation loop

%% Buoy centre and equation

cx = 1; % centre of gravity poisition in x-direction

cy = 0.035 + d; % centre of gravity position in y-direction

xc = r * cos(th) + cx ; % x values of buoy relative to poition in cartesian co-ordinates

yc = @(Y) r * sin(th) + cy + Y; % y values of buoy relative to poition in cartesian co-ordinates

centre_check = @(Y) cy + Y;

crc = @(Y) [xc;yc(Y)]; % Buoy geometry

top = @(Y) cy + r + Y; % highest point of buoy

bottom = @(Y) cy - r + Y; % lowest point of buoy

%% Wave surface elevation

xp = 0:inclambda:numwave*L; % horizontal distance [m] ** make a rounder number**

N =@(T) d + A*cos(k*(xp+cx) + w*T); % Surface of the wave (assuming Airy's linear wave theory)

wprofile = @(T) [xp ; N(T)]; % wave x and y co-ordinate

% water particle velocity

U1_s = @(T) ((pi*H)/Ti)*exp((2*pi*wl)/L)*cos(w*T); % water particle velocity in the surge (x-direction)

U2_h = @(T) ((pi*H)/Ti)*exp((2*pi*wl)/L)*sin(w*T); % water particle velocity in the heave (y-direction)

%% intersection of buoy and wave

%line

x_value(1:4) = cx;

dy = 0.2;

y_value = 0.2:dy:0.8;

line = [x_value;y_value];

select = @(Q,Z,S) Q(Z,S);

ptop = @(T) InterX(line,wprofile(T));

p11top = @(T) select(ptop(T), 1, 1);

p21top = @(T) select(ptop(T), 2, 1);

p =@(T,Y) InterX(crc(Y),wprofile(T)); % intersection points of wave and buoy done using the interX function

select = @(Q,Z,S) Q(Z,S);

tf = @(T,Y) isempty(p(T,Y)); % Is the matrix empty == true(1) or false (0) [ is there intersection points from the buoy and wave surface)

if tf(t(i),y(i)) == 1 && bottom(y(i)) > p21top(t(i)) %if intersection matrix is empty and bottom of buoy is above the wave surface

Fbuoy_h = @(T,Y)(T*Y)*0; % not in wave therefore no buoyancy

F_drag_h = @(T,Y,V) (T*Y*V)*0;

F_drag_h_wp = @(T,Y) (T*Y)*0;

Vol_sub = @(T,Y) (T*Y)*0;

F_Rd_h = @(T,Y)(T*Y)*0;

F_s_h = @(Y) S_k*(centre_check(Y)); % spring force in heave direction

Vol_percentage = @(T,Y) (T*Y)*0;

ifx(i) = 1

elseif tf(t(i),y(i)) == 1 && top(y(i)) < p21top(t(i))

A_sub =@(T,Y) pi*(r^2);

Vol_sub = @(T,Y) Vbuoy_Total;

Fbuoy_h = @(T,Y) (T*Y)*0 + Fbmax;

F_Rd_h = @(T,Y) (1/2)*Vol_sub(T,Y)*w*eh;

F_s_h = @(Y) S_k*(centre_check(Y)); % spring force in heave direction

F_drag_h = @(T,Y,V) 1/2*Cd*DSw*((V)*abs(V))*A_sub(T,Y);

F_drag_h_wp = @(T,Y) 1/2*Cd*DSw*(U2_h(T)*abs(U2_h(T)))*A_sub(T,Y);

Vol_percentage = @(T,Y) (T*Y)*0 +100;

ifx(i) = 2

elseif tf(t(i),y(i)) == 0 && centre_check(y(i)) > p21top(t(i)) %if intersection matrix is full and centre of buoy is above the wave surface

p11 = @(T,Y) select(p(T,Y), 1, 1);

p12 = @(T,Y) select(p(T,Y), 1, 2);

p21 = @(T,Y) select(p(T,Y), 2, 1);

p22 = @(T,Y) select(p(T,Y), 2, 2);

Distance_interX = @(T,Y) sqrt(((p12(T,Y) - p11(T,Y))^2) + ((p22(T,Y) - p21(T,Y))^2)); % distance between intersection point, Spherical cap diameter (2r)

Spherical_cap = @(T,Y)Distance_interX(T,Y)/2; % spherical cap base radius

Spherical_Distance =@(T,Y) sqrt((r^2) - (Spherical_cap(T,Y)^2)); % distance from centre of buoy to wave surface [m]

Cap_height =@(T,Y) r - Spherical_Distance(T,Y); % Spherical cap height relative to wave surface [m]

Vol_sub = @(T,Y)((pi*(Cap_height(T,Y)^2))/3)*((3*r)-Cap_height(T,Y)); % volume of spherical buoy submerged [m3]

Vol_percentage = @(T,Y)Vol_sub(T,Y)/Vbuoy_Total *100; % percentage of buoy submerged

A_sub =@(T,Y)pi*(Spherical_cap(T,Y)^2);

Fbuoy_h = @(T,Y) DSw*g*Vol_sub(T,Y); % buoyancy force expericenced by float [N]

F_drag_h = @(T,Y,V) 1/2*Cd*DSw*((V)*abs(V))*A_sub(T,Y);

F_Rd_h = @(T,Y) (1/2)*Vol_sub(T,Y)*w*eh;

F_s_h = @(Y) S_k*(centre_check(Y)); % spring force in heave direction

F_drag_h_wp = @(T,Y) 1/2*Cd*DSw*(U2_h(T)*abs(U2_h(T)))*A_sub(T,Y);

ifx(i) = 3

elseif tf(t(i),y(i)) == 0 && centre_check(y(i)) < p21top(t(i)) %if intersection matrix is full and centre of buoy is below the wave surface

p11 = @(T,Y) select(p(T,Y), 1, 1);

p12 = @(T,Y) select(p(T,Y), 1, 2);

p21 = @(T,Y) select(p(T,Y), 2, 1);

p22 = @(T,Y) select(p(T,Y), 2, 2);

Distance_interX = @(T,Y) sqrt(((p12(T,Y) - p11(T,Y))^2) + ((p22(T,Y) - p21(T,Y))^2)); % distance between intersection point, Spherical cap diameter (2r)

Spherical_cap = @(T,Y)Distance_interX(T,Y)/2; % spherical cap base radius

Spherical_Distance =@(T,Y) sqrt((r^2) - (Spherical_cap(T,Y)^2)); % distance from centre of buoy to wave surface [m]

Cap_height =@(T,Y) r - Spherical_Distance(T,Y); % Spherical cap height relative to wave surface [m]

Vol_sub = @(T,Y) ((Vbuoy_Total/2) + ((Vbuoy_Total/2) - ((pi*(Cap_height(T,Y)^2))/3)*((3*r)-Cap_height(T,Y)))); % volume of spherical buoy submerged [m3]

A_sub =@(T,Y)pi*(Spherical_cap(T,Y)^2);

Vol_percentage = @(T,Y)Vol_sub(T,Y)/Vbuoy_Total *100; % percentage of buoy submerged

Fbuoy_h = @(T,Y) DSw*g*Vol_sub(T,Y); % buoyancy force expericenced by float [N]

F_drag_h = @(T,Y,V) 1/2*Cd*DSw*((V)*abs(V))*A_sub(T,Y);

F_Rd_h = @(T,Y) (1/2)*Vol_sub(T,Y)*w*eh;

F_s_h = @(Y) S_k*(centre_check(Y)); % spring force in heave direction

F_drag_h_wp = @(T,Y) 1/2*Cd*DSw*(U2_h(T)*abs(U2_h(T)))*A_sub(T,Y);

ifx(i) = 4

end

added_mass_h = @(T,Y) DSw*Vol_sub(T,Y); %added mass for acceleration

Sum_forces_h = @(T,Y,V) Fbuoy_h(T,Y) - Weight + F_drag_h(T,Y,V);% + F_Rd_h(T,Y) + F_s_h(Y) + F_drag_h_wp(T,Y); % + added_mass_y(i);

%% Equation of motion to be intergrated using the Runge-Kutta 4th order method for Ordinary Diffferential Equation's (ODE's)

Accel_h = @(T,Y,V) ((Sum_forces_h(T,Y,V))/(abs(mass + added_mass_h(T,Y))));

k_1vh = Accel_h(t(i), y(i), v(i))*h;

k_1h = (v(i))*h;

k_2vh = Accel_h(t(i) + 0.5*h, y(i) + 0.5*k_1h, v(i) + 0.5*k_1vh)*h;

k_2h = (v(i) + (k_1vh/2))*h;

k_3vh = Accel_h((t(i) + 0.5*h),(y(i) + 0.5*k_2h) , v(i) + 0.5*k_2vh)*h;

k_3h = ((v(i)+(k_2vh/2))*h);

k_4vh = Accel_h((t(i)+h), (y(i)+k_3h), (v(i)+k_3vh))*h;

k_4h = (v(i) + k_3vh)*h;

n(i) = (1/6)*(k_1h + 2*k_2h + 2*k_3h + k_4h);

q(i) = (1/6)*(k_1vh + 2*k_2vh + 2*k_3vh + k_4vh);

y(i+1) = y(i) + n(i); % main equation cx

v(i+1) = v(i) + q(i); % main equation cy

%% Plots

vols(i) = Vol_percentage(t(i),y(i));

fsum(i) = Sum_forces_h(t(i),y(i),v(i));

tfx (i) = tf(t(i),y(i));

p21topp(i) = p21top(t(i));

bottomp(i) = bottom(y(i));

ip11(i) = p11(t(i), y(i));

ip12(i) = p12(t(i), y(i));

ip21(i) = p21(t(i), y(i));

ip22(i) = p22(t(i), y(i));

N = d + A*cos(k*(xp) + w*t(i)); % Surface of the wave (assuming Airy's linear wave theory)

xcn = r * cos(th) + cx; % x values of buoy relative to poition in cartesian co-ordinates

ycx =yc(y(i)); % y values of buoy relative to poition in cartesian co-ordinates

crcx = [xcn;ycx]; % Buoy geometry

figure(1)

title('motion of buoy on wave')

plot(xp,N,...

crcx(1,:),crcx(2,:))

ylim([0.25;0.7])

end

dpb
on 28 Apr 2019

Well, with all the code you did show, you didn't show us anything related to the error itself...namely the function that had the error.

However since the crystal ball does seem to have been repaired, and it even has a little bit of information to work with this time; let's put it to the test again-- :)

Index in position 2 exceeds array bounds.

Error in Floating_perfect_Save_Spare>@(Q,Z,S)Q(Z,S)

Error in Floating_perfect_Save_Spare>@(T,Y)select(p(T,Y),1,2)

Since the error message says you've exceeded an array bound and it's specifically the second argument to the function which is Z, one would presume the most likely cause is you have a loop over the size of the array that is looking at positions Z(i),Z(i+1) and when i reaches the size of the array Z, then Z(i+1) is out of bounds.

Very common logic error to make; simply make the upper loop bound one less than the array size in the particular dimension.

Of course, you need to also think about your intersection search logic -- if it is at some boundary, can you actually reach the intersection itself by the search you are using? Work out with pencil and paper what you're actually looking at and see for the end conditions...

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.