Can you help me write this question in finite volume method? This is my first time writing FVM. I wrote using FDM below :(

7 views (last 30 days)
%Upwind scheme
n=320;
X=n+1;
Y=n+1;
u=0.5;
v=0.5;
Lx=4;
Ly=4;
dX=Lx/n;
dY=Ly/n;
c(:,:,1)=zeros(Y,X);
c(:,:,2)=zeros(Y,X);
for a=(2+0.9/dX):1:(1.1/dX)
for b=(2+0.9/dY):1:(1.1/dY)
x=(a-1)*dX;
y=(b-1)*dY;
c(a,b,1)=sin(5*pi*(x-0.9))*sin(5*pi*(y-0.9));
end
end
dT=0.001;
for T=1:1:4
for n=2:1:(round(T/dT)+1)
for j=2:1:(Lx/dX)
for i=2:1:(Ly/dY)
c(i,j,2)=c(i,j,1)+(v*dT/dX)*(c(i-1,j,1)-c(i,j,1))+(u*dT/dY)*(c(i,j-1,1)-c(i,j,1));
end
end
c(:,:,1)=c(:,:,2);
end
dT=dT+0.001;
fig = figure();
mesh(0:dX:Lx,0:dY:Ly,c(:,:,2));
xlabel('X');
ylabel('Y');
zlabel('c');
grid on;
axis([0 4 0 4 -1 1]);
end

Answers (1)

Brahmadev
Brahmadev on 9 Feb 2024
The MATLAB code for First order upwind seems to be on the right track, here are some more things to consider:
  1. The time loop has an increment dT=dT+0.001; at the end, which doesn't make sense in this context since dT is your time step and should remain constant as defined at the beginning.
  2. The boundary conditions are not enforced in your code. You should set the edges of the computational domain to zero at each time step.
Here is a revised version of the code that takes care of these considerations:
% Define the domain and grid size
n = 320;
Lx = 4;
Ly = 4;
dx = Lx / n;
dy = Ly / n;
x = linspace(0, Lx, n+1);
y = linspace(0, Ly, n+1);
[X, Y] = meshgrid(x, y);
% Define the time domain
dt = 0.001;
t_final = 4;
Nt = round(t_final / dt);
% Define the initial condition
c = zeros(n+1, n+1);
for i = 2:n
for j = 2:n
if (x(i) > 0.9 && x(i) < 1.1) && (y(j) > 0.9 && y(j) < 1.1)
c(j, i) = sin(5 * (x(i) - 0.9)) * sin(5 * (y(j) - 0.9));
end
end
end
% Define velocities
u = 0.5;
v = 0.5;
% FVM implementation
for t = 1:Nt
c_new = c;
% Compute fluxes and update the scalar field c
for i = 2:n
for j = 2:n
% Compute fluxes using upwind scheme
flux_x = u * (c(j, i) - c(j, i-1)) / dx;
flux_y = v * (c(j, i) - c(j-1, i)) / dy;
% Update c based on net fluxes
c_new(j, i) = c(j, i) - dt * (flux_x + flux_y);
end
end
% Apply boundary conditions
c_new(1, :) = 0; % Bottom boundary
c_new(end, :) = 0; % Top boundary
c_new(:, 1) = 0; % Left boundary
c_new(:, end) = 0; % Right boundary
c = c_new;
% Plot at specific time steps
if mod(t, Nt / 4) == 0 || t == 1
figure;
surf(X, Y, c, 'EdgeColor', 'none');
title(sprintf('FVM with Upwind Scheme at t = %.3f', t * dt));
xlabel('x');
ylabel('y');
zlabel('c');
view(3);
colorbar;
pause(0.1); % To visualize the evolution
end
end
Hope this helps in resolving your query!

Categories

Find more on Computational Fluid Dynamics (CFD) in Help Center and File Exchange

Products

Community Treasure Hunt

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

Start Hunting!