Simulation of Faraday Waves with Matlab

Dear all,

I have a matlab code described here below which is running but which is incomplete. I would like to visualize the surface displacement of the fluid as a function of the x and y coordinates, but I get only a flat surface. Any help is welcomed !

% Define simulation parameters

rho = 1000; % fluid density [kg/m^3]

mu = 1e-3; % fluid viscosity [Pa*s]

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]

% Set up grid

Nx = 100;

Ny = 50;

x = linspace(0, L, Nx);

y = linspace(0, W, Ny);

[X, Y] = meshgrid(x, y);

% Define simulation parameters

rho = 1000; % fluid density [kg/m^3]

mu = 1e-3; % fluid viscosity [Pa*s]

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]

% Set up grid

Nx = 100;

Ny = 50;

x = linspace(0, L, Nx);

y = linspace(0, W, Ny);

[X, Y] = meshgrid(x, y);

% 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]

% Run simulation

for i = 2:length(t)

% Compute acceleration at current time step

[a, b] = acceleration(u, v, eta, rho, mu, g, A, omega, t(i));

% Update velocity and displacement using Euler method

u = u + a*dt;

v = v + b*dt;

eta = eta + v*dt;

end

% Visualize results

figure;

surf(X, Y, eta);

xlabel('x');

ylabel('y');

zlabel('displacement');

title('Faraday waves');

% Define function to compute acceleration

function [a, b] = acceleration(u, v, eta, rho, mu, g, A, omega, t)

% Compute acceleration using Navier-Stokes equations

a = -rho*g*eta - rho*A*omega^2*sin(omega*t);

b = -rho*g*eta - rho*A*omega^2*cos(omega*t);

end

Answer 1

[Sulaymon Eshkabilov]

Here is a corrected code of your exercise:

clearvars; clc
% Define simulation parameters
rho = 1000;    % fluid density [kg/m^3]
mu = 1e-3;     % fluid viscosity [Pa*s]
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]
% Set up grid
Nx = 100;
Ny = 50;
% Set up time-stepping
dt = 0.005;    % time step [s]
tmax = 10;     % maximum time [s]
t = 0:dt:(dt*(Nx-1)); % 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]
% Run simulation
 for ii = 2:Nx
    for jj = 2:Ny
  % Compute acceleration at current time step
  a = -rho*g*eta(ii-1, jj-1) - rho*A*omega^2*sin(omega*t(ii));
  b = -rho*g*eta(ii-1, jj-1) - rho*A*omega^2*cos(omega*t(ii));
  % Update velocity and displacement using Euler method
  u(ii,jj) = u(ii-1, jj-1) + a*dt;
  v(ii,jj) = v(ii-1, jj-1) + b*dt;
  eta(ii,jj) = eta(ii-1,jj-1) + v(ii, jj)*dt;
    end
end
% Visualize results
figure;
surf(u, v, eta);
xlabel('u');
ylabel('v');
zlabel('displacement');
title('Faraday waves');