New code to see why it doesn't work

Question

Thierry 2023-1-12

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1893100-new-code-to-see-why-it-doesn-t-work

编辑： Khalid 2023-1-18

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

6 个评论
显示 4更早的评论隐藏 4更早的评论

Thierry 2023-1-12

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

Thierry 2023-1-13

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 2023-1-13

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 2023-1-16

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1893100-new-code-to-see-why-it-doesn-t-work#answer_1149470

编辑：Khalid 2023-1-18

在 MATLAB Online 中打开

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

New code to see why it doesn't work

6 个评论
显示 4更早的评论隐藏 4更早的评论

回答（1 个）

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

New code to see why it doesn't work

6 个评论 显示 4更早的评论隐藏 4更早的评论

回答（1 个）

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

6 个评论
显示 4更早的评论隐藏 4更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论