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

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

Answer 1

[VBBV]

% 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

[VBBV]

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]

thank you very much VBBV, it now works!