New code to see why it doesn't work

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);
% Compute pressure gradient force
[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

[Torsten]

"it doesn't work properly" is a very vague statement.

A litlle bit of explanation about what your code does and where you encounter difficulties would be helpful.

[Steven Lord]

What does "doesn't work properly" mean in this context?

• Do you receive warning and/or error messages? If so the full and exact text of those messages (all the text displayed in orange and/or red in the Command Window) may be useful in determining what's going on and how to avoid the warning and/or error.
• Does it do something different than what you expected? If so, what did it do and what did you expect it to do?
• Did MATLAB crash? If so please send the crash log file (with a description of what you were running or doing in MATLAB when the crash occured) to Technical Support so we can investigate.

[Thierry]

Dear Torsten, dear Steven,

Sorry for my laconic explanation about my difficulties I encountered with this code. Here some useful explanations about the code and what I want to achieve:

this code should compute the movement of several floating particles on the surface of a fluid (water for instance) subjected to a parametric vertical oscillation, the so called Faraday waves. The task is known to be difficult because the wave flow is nonlinear and can be described by the Navier-Stokes equations (fundamental equations used in hydrodynamics to describe the movement of fluids (laminar and turbulent)). Instead of it, I use here three forces: 1) the buoyancy force, 2) the drag force, and 3) the pressure gradient force. From the sum of these three forces, I get the acceleration acting on each floater (j) at each time step (i). The goal of the program it to see the evolution of the system depending of time and space.

When I run the program, I get the following message in red:

"Arrays have incompatible sizes for this operation.

Error in GoOn2>acceleration (line 85)

F_total = F_b + F_d + F_p;

Error in GoOn2 (line 41)

[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));

Related documentation"

I expect to visualize the plot of the computation with different snapshots at different times, with the displacement (in z direction) of the different floaters at different time. I get no figure at all. MATLAB doesn't crash.

I hope these explanations are sufficient.

Thank you again!

Thierry

[Torsten]

What shall we say ?

Make the sizes of F_b, F_d and F_p equal so that they can be added. You see the sizes in your code above - I printed them for you.

[Thierry]

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]

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);

% Compute pressure gradient force

[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

Answer 1

[Khalid]

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.