% 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
