Simulation of Faraday Waves with Matlab

15 次查看（过去 30 天）
Thierry
Thierry 2023-1-7
编辑： Sulaymon Eshkabilov 2023-1-7
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

请先登录，再回答此问题。

回答（1 个）

Sulaymon Eshkabilov
编辑：Sulaymon Eshkabilov 2023-1-7
Ran in:
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');

类别

在 Help Center 和 File Exchange 中查找有关 Fluid Dynamics 的更多信息

标签

产品


版本

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by