# New code to see why it doesn't work

3 views (last 30 days)
Thierry on 12 Jan 2023
Edited: Khalid on 18 Jan 2023
Here is a new code I created, but unfortunately it doesn't work properly. Could anyone help me the program to work? Thank you a lot!
% Define simulation parameters
rho = 1000; % fluid density [kg/m^3]
mu = 1e-3; % fluid viscosity [Pas]
g = 9.81; % acceleration due to gravity [m/s^2]
L = 1; % length of container [m]
W = 0.5; % width of container [m]
A = 0.1; % amplitude of oscillation [m]
omega = 2*pi; % frequency of oscillation [rad/s]
num_floaters = 10; % number of floaters
rho_f = 500*ones(1,num_floaters); % floater density [kg/m^3]
V = 0.01*ones(1,num_floaters); % floater volume [m^3]
% Set up grid
Nx = 100;
Ny = 50;
x = linspace(0, L, Nx); % intervalle 0 - L divisé en Nx segments
y = linspace(0, W, Ny);
[X, Y] = meshgrid(x, y); %creates two matrices X and Y of the same size that represent a 2-dimensional grid. The %matrix X has the same value as x in each row, and the matrix Y has the same value as y in each column. In this %specific case, the x variable is a 1-dimensional array of Nx evenly spaced values between 0 and L, and the y %variable is a 1-dimensional array of Ny evenly spaced values between 0 and W. The purpose of creating X and Y %matrices is to have a set of coordinates on a 2-dimensional grid with the x-coordinates ranging from 0 to L and y-%coordinates ranging from 0 to W, and this will be used later in the code for visualization and computations.
dx = L/Nx; % grid spacing in x direction
dy = W/Ny; % grid spacing in y direction
% Set up time-stepping
dt = 0.001; % time step [s]
tmax = 10; % maximum time [s]
t = 0:dt:tmax; % time vector
% Set up initial conditions
u = zeros(Ny, Nx); % initial x velocity [m/s]
v = zeros(Ny, Nx); % initial y velocity [m/s]
eta = zeros(Ny, Nx); % initial displacement [m]
x_floater = L*rand(1,num_floaters); % initial x position of floater [m]
y_floater = W*rand(1,num_floaters); % initial y position of floater [m]
a_floater = zeros(size(X));
b_floater = zeros(size(Y));
% Run simulation
for i = 2:length(t)
% Compute acceleration at current time step
for j = 1:num_floaters
[a(j), b(j)] = acceleration(u, v, dx, dy, X, Y, eta, rho, mu, g, A, omega, t(i), rho_f(j), V(j), x_floater(j), y_floater(j));
end
a_floater = sum(a);
b_floater = sum(b);
% Update velocity and displacement using Euler method
u = u + a_floater*dt;
v = v + b_floater*dt;
eta = eta + u*dt + v*dt;
% Update position of floater
x_floater = x_floater + u(round(y_floater/dy), round(x_floater/dx))*dt;
y_floater = y_floater + v(round(y_floater/dy), round(x_floater/dx))*dt;
% Check if floater is within bounds of container
x_floater(x_floater < 0) = 0;
x_floater(x_floater > L) = L;
y_floater(y_floater < 0) = 0;
y_floater(y_floater > W) = W;
% Plot displacement and floater position
figure(1);
surf(X, Y, eta);
hold on;
plot3(x_floater, y_floater, eta(round(y_floater/dy), round(x_floater/dx)), 'r.', 'MarkerSize', 15);
hold off;
xlabel('x [m]');
ylabel('y [m]');
zlabel('displacement [m]');
title(['time = ', num2str(t(i)), 's']);
drawnow;
end
ans = 1×2
1 1
ans = 1×2
50 100
ans = 1×2
1 2
Arrays have incompatible sizes for this operation.

Error in solution>acceleration (line 82)
F_total = F_b + F_d + F_p;
% Define function to compute acceleration
function [a, b] = acceleration(u, v, dx, dy, X, Y, eta, rho, mu, g, A, omega, t, rho_f, V, x_floater, y_floater)
% Compute gradient of velocity field
[du_dx, du_dy] = gradient(u, dx, dy);
[dv_dx, dv_dy] = gradient(v, dx, dy);
% Compute pressure field
p = -rho*g*eta;
% Compute buoyancy force
F_b = rho_f*g*V;
% Compute drag force
F_d = -0.5*rho_f*V*(du_dx + dv_dy);
[dp_dx, dp_dy] = gradient(p, dx, dy);
F_p = [dp_dx(round(y_floater/dy), round(x_floater/dx)), dp_dy(round(y_floater/dy), round(x_floater/dx))];
% Compute total force on floater
size(F_b)
size(F_d)
size(F_p)
F_total = F_b + F_d + F_p;
% Compute acceleration
a = F_total(1)/(rho_f*V);
b = F_total(2)/(rho_f*V);
b = b - A*omega*sin(omega*t);
end
Thierry on 13 Jan 2023
Hello, thank you for your help. Since I did not develop this code by myself but used artificial intelligence instead, I understand the meaning of the code but I am still new to Matlab computing.
So, what should I add explicitly in the code in order to make the sizes of F_b, F_d and F_p equal so that they can be added? I tried with the following lines of code inserted just before the calculation of F_total:
F_b = zeros(size(X));
F_d = zeros(size(X));
F_p = zeros(size(X));:
But then I get the following message:
Arrays have incompatible sizes for this operation.
Error in GoOn2 (line 57)
y_floater = y_floater + v(round(y_floater/dy), round(x_floater/dx))*dt;
How can I get rid of this error message?
Thank you again for your valuable help!
Thierry
Thierry on 13 Jan 2023
Here is the last version of the code:
% Define simulation parameters
rho = 1000; % fluid density [kg/m^3]
mu = 1e-3; % fluid viscosity [Pas]
g = 9.81; % acceleration due to gravity [m/s^2]
L = 1; % length of container [m]
W = 0.5; % width of container [m]
A = 0.1; % amplitude of oscillation [m]
omega = 2*pi; % frequency of oscillation [rad/s]
num_floaters = 10; % number of floaters
rho_f = 500*ones(1,num_floaters); % floater density [kg/m^3]
V = 0.01*ones(1,num_floaters); % floater volume [m^3]
% Set up grid
Nx = 100;
Ny = 50;
x = linspace(0, L, Nx); % intervalle 0 - L divisé en Nx segments
y = linspace(0, W, Ny);
[X, Y] = meshgrid(x, y); %creates two matrices X and Y of the same size that represent a 2-dimensional grid. The matrix X has the same value as x in each row, and the matrix Y has the same value as y in each column. In this specific case, the x variable is a 1-dimensional array of Nx evenly spaced values between 0 and L, and the y variable is a 1-dimensional array of Ny evenly spaced values between 0 and W. The purpose of creating X and Y matrices is to have a set of coordinates on a 2-dimensional grid with the x-coordinates ranging from 0 to L and y-coordinates ranging from 0 to W, and this will be used later in the code for visualization and computations.
dx = L/Nx; % grid spacing in x direction
dy = W/Ny; % grid spacing in y direction
% Set up time-stepping
dt = 0.001; % time step [s]
tmax = 10; % maximum time [s]
t = 0:dt:tmax; % time vector
% Set up initial conditions
u = zeros(Ny, Nx); % initial x velocity [m/s]
v = zeros(Ny, Nx); % initial y velocity [m/s]
eta = zeros(Ny, Nx); % initial displacement [m]
x_floater = L*rand(1,num_floaters); % initial x position of floater [m]
y_floater = W*rand(1,num_floaters); % initial y position of floater [m]
a_floater = zeros(size(X));
b_floater = zeros(size(Y));
% Run simulation
for i = 2:length(t)
% Compute acceleration at current time step
for j = 1:num_floaters
[a(j), b(j)] = acceleration(u, v, dx, dy, X, Y, eta, rho, mu, g, A, omega, t(i), rho_f(j), V(j), x_floater(j), y_floater(j));
end
a_floater = sum(a);
b_floater = sum(b);
% Update velocity and displacement using Euler method
u = u + a_floater*dt;
v = v + b_floater*dt;
eta = eta + u*dt + v*dt;
% Update position of floater
% % % x_floater = zeros(size(Y));
% % % y_floater = zeros(size(X));
% % % u(round(y_floater/dy), round(x_floater/dx)) = zeros(size(X));
% % % v(round(y_floater/dy), round(x_floater/dx)) = zeros(size(X));
x_floater = x_floater + u(round(y_floater/dy), round(x_floater/dx))*dt;
y_floater = y_floater + v(round(y_floater/dy), round(x_floater/dx))*dt;
% Check if floater is within bounds of container
x_floater(x_floater < 0) = 0;
x_floater(x_floater > L) = L;
y_floater(y_floater < 0) = 0;
y_floater(y_floater > W) = W;
% Plot displacement and floater position
figure(1);
surf(X, Y, eta);
hold on;
plot3(x_floater, y_floater, eta(round(y_floater/dy), round(x_floater/dx)), 'r.', 'MarkerSize', 15);
hold off;
xlabel('x [m]');
ylabel('y [m]');
zlabel('displacement [m]');
title(['time = ', num2str(t(i)), 's']);
drawnow;
end
% Define function to compute acceleration
function [a, b] = acceleration(u, v, dx, dy, X, Y, eta, rho, mu, g, A, omega, t, rho_f, V, x_floater, y_floater)
% Compute gradient of velocity field
[du_dx, du_dy] = gradient(u, dx, dy);
[dv_dx, dv_dy] = gradient(v, dx, dy);
% Compute pressure field
p = -rho*g*eta;
% Compute buoyancy force
F_b = rho_f*g*V;
% Compute drag force
F_d = -0.5*rho_f*V*(du_dx + dv_dy);
[dp_dx, dp_dy] = gradient(p, dx, dy);
F_p = [dp_dx(round(y_floater/dy), round(x_floater/dx)), dp_dy(round(y_floater/dy), round(x_floater/dx))];
% Compute total force on floater
F_b = zeros(size(X));
F_d = zeros(size(X));
F_p = zeros(size(X));
F_total = F_b + F_d + F_p;
% Compute acceleration
a = F_total(1)/(rho_f*V);
b = F_total(2)/(rho_f*V);
b = b - A*omega*sin(omega*t);
end

Khalid on 16 Jan 2023
Edited: Khalid on 18 Jan 2023
Hello,
It is my understanding that you encountered the error “Arrays have incompatible sizes for this operation” while executing the code.
After line 56 gets executed, 'x_floater' becomes a matrix of size 10 x 10 which was of size 1 x 10 before execution of line 56. Since 'x_floater's size is 10x10 in line 57, you are receiving the error saying incompatible array sizes.
As a simple workaround for this issue:
x_new_floater = x_floater + u(round(y_floater/dy), round(x_floater/dx))*dt;
y_new_floater = y_floater + v(round(y_floater/dy), round(x_floater/dx))*dt;
x_floater = x_new_floater;
y_floater = y_new_floater;
Here, we are renaming the updated floaters 'x_floater' and ‘y_floater' to 'x_new_floater' and 'y_new_floater' and using the old floater values to compute these updates, making sure that the dimensions of 'x_floater' and 'y_floater' remain same.

### Categories

Find more on Fluid Mechanics in Help Center and File Exchange

R2022b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!