Clear Filters
Clear Filters

Simulation of Faraday Waves with Matlab

15 views (last 30 days)
Thierry
Thierry on 7 Jan 2023
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

Answers (1)

Sulaymon Eshkabilov
Sulaymon Eshkabilov on 7 Jan 2023
Edited: Sulaymon Eshkabilov on 7 Jan 2023
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');

Categories

Find more on Fluid Dynamics in Help Center and File Exchange

Products


Release

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!