New code to see why it doesn't work
2 次查看(过去 30 天)
显示 更早的评论
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
% 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
回答(1 个)
Khalid
2023-1-16
编辑:Khalid
2023-1-18
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 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Fluid Mechanics 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!