# My simulation doesn't take into account the different height of the Faraday waves?

4 views (last 30 days)
Thierry on 11 Dec 2022
Commented: Thierry on 8 Jan 2023
Could you please tell me why the following plot doesn't give a correct plotting of the different heights of the Faraday waves? Thanks a lot!
% Define physical parameters
g = 9.81; % acceleration due to gravity
d = 0.1; % fluid depth
omega = 0.1; % forcing frequency
% Define grid parameters
Nx = 100; % number of x grid points
Ny = 100; % number of y grid points
Lx = 1; % length of x dimension
Ly = 1; % length of y dimension
dx = Lx/Nx; % grid spacing in x direction
dy = Ly/Ny; % grid spacing in y direction
[x, y] = meshgrid(dx/2:dx:Lx-dx/2, dy/2:dy:Ly-dy/2); % grid
% Define initial conditions
u = zeros(Nx, Ny); % x-velocity
v = zeros(Nx, Ny); % y-velocity
h = zeros(Nx, Ny); % surface height
% Define forcing function
F = @(t) 0.1*sin(omega*t);
% Define time-stepping parameters
dt = 0.001; % time step
tmax = 10; % maximum time
t = 0:dt:tmax; % time vector
% Perform time-stepping
for ii = 1:length(t)
% Loop through grid points, updating surface height and velocities
for jj = 1:Nx
for kk = 1:Ny
% Calculate acceleration using equations of motion
a_x = -g*h(jj,kk)*sin(h(jj,kk)) - g/2*(u(jj,kk)^2 + v(jj,kk)^2) + F(t(ii));
a_y = -g*h(jj,kk)*sin(h(jj,kk)) - g/2*(u(jj,kk)^2 + v(jj,kk)^2) + F(t(ii));
% Update velocities
u(jj,kk) = u(jj,kk) + a_x*dt;
v(jj,kk) = v(jj,kk) + a_y*dt;
% Update surface height
h(jj,kk) = h(jj,kk) + dt*(u(jj,kk)*v(jj,kk) - g/2*h(jj,kk)^2);
end
end
end
% Plot results
figure;
mesh(x, y, h);
xlabel('x');
ylabel('y');
zlabel('surface height');
##### 0 CommentsShow -2 older commentsHide -2 older comments

Sign in to comment.

### Accepted Answer

VBBV on 7 Jan 2023
% Define physical parameters
g = 9.81; % acceleration due to gravity
d = 0.1; % fluid depth
omega = 0.1; % forcing frequency
% Define grid parameters
Nx = 100; % number of x grid points
Ny = 100; % number of y grid points
Lx = 1; % length of x dimension
Ly = 1; % length of y dimension
dx = Lx/Nx; % grid spacing in x direction
dy = Ly/Ny; % grid spacing in y direction
[x, y] = meshgrid(dx/2:dx:Lx-dx/2, dy/2:dy:Ly-dy/2); % grid
% Define initial conditions
u = zeros(Nx, Ny); % x-velocity
v = zeros(Nx, Ny); % y-velocity
h = zeros(Nx, Ny); % surface height
% Define forcing function
F = @(t) 0.1*sin(omega*t);
% Define time-stepping parameters
dt = 0.1; % time step
tmax = 1; % maximum time
t = 0:dt:tmax; % time vector
% Perform time-stepping
for ii = 1:length(t)
% Loop through grid points, updating surface height and velocities
for jj = 1:Nx-1
for kk = 1:Ny-1
% Calculate acceleration using equations of motion
a_x = -g*h(jj,kk)*sin(h(jj,kk)) - g/2*(u(jj,kk)^2 + v(jj,kk)^2) + F(t(ii));
a_y = -g*h(jj,kk)*sin(h(jj,kk)) - g/2*(u(jj,kk)^2 + v(jj,kk)^2) + F(t(ii));
% Update velocities
u(jj+1,kk+1) = u(jj,kk) + a_x*dt;
v(jj+1,kk+1) = v(jj,kk) + a_y*dt;
% Update surface height
h(jj+1,kk+1) = h(jj,kk) + dt*(u(jj,kk)*v(jj,kk) - g/2*h(jj,kk)^2);
end
end
nexttile
mesh(x, y, h);
xlabel('x');
ylabel('y');
zlabel('h');
title(sprintf('t=%0.1ds',t(ii))); % at different times
%view(45,90)
end
##### 2 CommentsShow NoneHide None
VBBV on 7 Jan 2023
Edited: VBBV on 7 Jan 2023
update the velocities at each dependent grid step, as
u(jj+1,kk+1) = u(jj,kk) + a_x*dt;
v(jj+1,kk+1) = v(jj,kk) + a_y*dt;
% Update surface height
h(jj+1,kk+1) = h(jj,kk) + dt*(u(jj,kk)*v(jj,kk) - g/2*h(jj,kk)^2);
it seems you were trying to compute wave heights at discrete time steps independently
Thierry on 8 Jan 2023
thank you very much VBBV, it now works!

Sign in to comment.

### Categories

Find more on Thermal Analysis 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!