CFD Online Discussion Forums - Matlab code for Finite Volume Method in 2D

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- Main CFD Forum (https://www.cfd-online.com/Forums/main/)

- - Matlab code for Finite Volume Method in 2D (https://www.cfd-online.com/Forums/main/144808-matlab-code-finite-volume-method-2d.html)

Matlab code for Finite Volume Method in 2D

Dear Forum members,

I recently begun to learn about basic Finite Volume method, and I am trying to apply the method to solve the following 2D continuity equation on the cartesian grid

$\rho_t + \nabla \cdot (\rho \textbf{u}) = 0$ with initial condition $\rho(x,y,0) = 1$

For simplicity and interest, I take $\textbf{u} = -\nabla \phi$ , where $\phi$ is the distance function given by $\phi = \sqrt({(x-\frac{1}{2})^2+ (y-\frac{1}{2})^2})$ so that all the density is concentrated near the point $(\frac{1}{2},\frac{1}{2})$ after sufficiently long enough time.

To solve this problem using the Finite Volume Method, I have written the matlab code (with uniform grid in x and y).

Code:

% Create Grid and number of cells

a = 0; 

b = 1; 

N = 10;

% Define edges

x_edges = linspace(a,b,N+1);

y_edges = linspace(a,b,N+1);

%Define distance between edges

delta_x = x_edges(2) - x_edges(1);

delta_y = y_edges(2) - y_edges(1);

%Define cell centers

x_centers = a+delta_x/2 : delta_x : b;

y_centers = a+delta_y/2 : delta_y : b;

[X,Y] = meshgrid(x_centers,y_centers);      % 2d arrays of x,y center values

X = X';                     % transpose so that X(i,j),Y(i,j) are

Y = Y';                     % (i,j) cell.

%Initialize solution array

U = zeros(N,N);

 

% Fill the solution array with initial data

U = X*0+1;

%define distance function

phi = sqrt((X-1/2).^2+(Y-1/2).^2);

gradphix = (X-1/2)./phi;

gradphiy = (Y-1/2)./phi;

% %reverse sign so that the gradient points inwards

v1 = -gradphix;

v2 = -gradphiy;

 

%Define timestep

delta_t = 0.001;

F = v1.*U;

G = v2.*U;

%Store initial data for later use

ic = U;

timesteps = 20;

for k = 1:timesteps

for j = 1:N % denote row position

 

     %denote column position

     for i = 1:N

 

        %Top row

        if (j ==N) 

            %Down flux

            G_down =  (1/delta_y)*(G(i,j)-G(i,j-1));

 

            %Nothing is coming down from above

            G_up =  0;

 

        %Bottom row    

        elseif (j == 1)

 

            %Nothing is coming in from below

            G_down = 0;

 

            %Up flux

            G_up = (1/delta_y)*(G(i,j+1)-G(i,j));

 

        %Rows in the middle

        else

 

            G_up = (1/delta_y)*(G(i,j+1)-G(i,j));

            G_down =  (1/delta_y)*(G(i,j)-G(i,j-1));

 

        end

 

 

        %Sideway fluxes

            if i == 1 

                F_left = 0;

                F_right = (1/delta_x)*(F(i+1,j)-F(i,j));

            elseif i == N

                F_left = (1/delta_x)*(F(i,j)-F(i-1,j));

                F_right = 0;

            else

                F_right = (1/delta_x)*(F(i+1,j)-F(i,j));

                F_left = (1/delta_x)*(F(i,j)-F(i-1,j));

            end

 

        Unew(i,j) = U(i,j) - (delta_t)*(F_right-F_left)-(delta_t)*(G_up - ...

             G_down);

     end

end

U = Unew;

surf(X,Y,Unew)

pause(0.1);

end

In my code, I have tried to implement a fully discrete flux-differencing method as on pg 440 of Randall LeVeque's Book "Finite Volume Methods for Hyperbolic Problems".

My code does not do its job, and I believe that there is something wrong with how I calculate my Fluxes through the four sides of my rectangular cell. ( F_right,F_left,G_up,G_down), and how I update my cell average after $\Delta t$ .

Could anyone give me some help?
I have spent the last two weeks working on this problem, but I am still stuck and do not know how to fix it
Also, is there a way to get rid of for-loops in my code?

My apologies if my question is too elementary for this forum!

Hello coagmento,

I have been working on developing a k-epsilon tubulence model in MatLab. Though my model has not been yet completed, I would like to know more more about your works.

Can we talk more about this topic? :)

My email address is : intishar.rahman@gmail.com

% From Michael
% 1)Your main error is the flux ! see below and compare
% 2) Secondly use reflection b.c.'s.
% 3) The scheme below is actually FTCS which is unstable for convection
% alone !!! so beware...
% Let me know ! what you is your project ?
% Create Grid and number of cells
a = 0;
b = 1;
N = 40;
% Define edges
x_edges = linspace(a,b,N+1);
y_edges = linspace(a,b,N+1);
%Define distance between edges
delta_x = x_edges(2) - x_edges(1);
delta_y = y_edges(2) - y_edges(1);
%Define cell centers
x_centers = a+delta_x/2 : delta_x : b;
y_centers = a+delta_y/2 : delta_y : b;
[X,Y] = meshgrid(x_centers,y_centers); % 2d arrays of x,y center values
X = X'; % transpose so that X(i,j),Y(i,j) are
Y = Y'; % (i,j) cell.
%Initialize solution array
U = zeros(N,N); Unew = zeros(N,N);

% Fill the solution array with initial data
U = X*0+1;
%define distance function
phi = sqrt((X-1/2).^2+(Y-1/2).^2);
gradphix = (X-1/2)./phi;
gradphiy = (Y-1/2)./phi;
% %reverse sign so that the gradient points inwards
v1 = -gradphix;
v2 = -gradphiy;

%Define timestep
delta_t = 0.001;
F = v1.*U;
G = v2.*U;
%Store initial data for later use
ic = U;
timesteps = 20;
for k = 1:timesteps
for j = 1:N % denote row position

%denote column position
for i = 1:N

%Top row
if (j ==N)
%Down flux
G_down = (1/delta_y)*0.5*(G(i,j)+G(i,j-1));

%Nothing is coming down from above
% G_up = 0;
G_up = G_down;

%Bottom row
elseif (j == 1)

%Nothing is coming in from below
% G_down = 0;

%Up flux
G_up = (1/delta_y)*0.5*(G(i,j+1)+G(i,j));
G_down =G_up;
%Rows in the middle
else

G_up = (1/delta_y)*0.5*(G(i,j+1)+G(i,j));
G_down = (1/delta_y)*0.5*(G(i,j)+G(i,j-1));

end

%Sideway fluxes
if i == 1
% F_left = 0;
F_right = (1/delta_x)*0.5*(F(i+1,j)+F(i,j));
F_left =F_right;
elseif i == N
F_left = (1/delta_x)*0.5*(F(i,j)+F(i-1,j));
% F_right = 0;
F_right =F_left;
else
F_right = (1/delta_x)*0.5*(F(i+1,j)+F(i,j));
F_left = (1/delta_x)*0.5*(F(i,j)+F(i-1,j));
end

Unew(i,j) = U(i,j) - (delta_t)*(F_right-F_left)-(delta_t)*(G_up - ...
G_down);
end
end
U = Unew;
surf(X,Y,Unew)
pause(0.1);
end

matlab code for Roe lineaarization

I am writing the code of Roe linearization for Shallow water equation.
I don't know how to implement it in matlab.
Is it possible to help me with that?
Thank you.

I am trying to write a Matlab program for a 1D unsteady conduction equation(parabolic pde).

i used finite volume method explicit scheme to solve the problem for a time=t+1.

I need idea on a Matlab code that would provide future iterations.

i am new to this field.

\thanks.

Finite volume 2d coaxial ground heat exchanger

Hi everyone
I am phd student, and i am stuck in bulding the program for simulation of the coaxial ground heat exchanger,i need some help in defining the programme skull,i will appreciate any help thank you all.