CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > General Forums > Main CFD Forum

Matlab code for Finite Volume Method in 2D

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree5Likes
  • 3 Post By coagmento
  • 1 Post By mgedwards
  • 1 Post By AliTalbi

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   November 22, 2014, 01:32
Default Matlab code for Finite Volume Method in 2D
  #1
New Member
 
David
Join Date: Nov 2014
Posts: 2
Rep Power: 0
coagmento is on a distinguished road
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 [0,1] x [0,1]

\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!
ramakant, mend4x and nasirhossain like this.
coagmento is offline   Reply With Quote

Old   December 15, 2014, 04:09
Thumbs up
  #2
New Member
 
Md. Intishar's Avatar
 
Md. Intishar Rahman
Join Date: Nov 2014
Location: Bangladesh
Posts: 6
Rep Power: 11
Md. Intishar is on a distinguished road
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
Md. Intishar is offline   Reply With Quote

Old   August 27, 2015, 05:45
Smile Corrected code
  #3
New Member
 
mgedwards
Join Date: Aug 2015
Posts: 2
Rep Power: 0
mgedwards is on a distinguished road
% 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
Alisha likes this.
mgedwards is offline   Reply With Quote

Old   March 21, 2016, 17:30
Default matlab code for Roe lineaarization
  #4
New Member
 
Join Date: Mar 2016
Posts: 2
Rep Power: 0
sisy is on a distinguished road
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.
sisy is offline   Reply With Quote

Old   March 26, 2016, 04:44
Default hello
  #5
New Member
 
leonard Isaac
Join Date: Mar 2016
Posts: 1
Rep Power: 0
leoimewore is on a distinguished road
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.
leoimewore is offline   Reply With Quote

Old   October 7, 2022, 13:47
Default Finite volume 2d coaxial ground heat exchanger
  #6
New Member
 
Okay
Join Date: Oct 2022
Posts: 1
Rep Power: 0
AliTalbi is on a distinguished road
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.
nasirhossain likes this.
AliTalbi is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
my new matlab code houkensjtu Main CFD Forum 3 July 1, 2016 16:28
alphaEqn.H in twoPhaseEulerFoam cheng1988sjtu OpenFOAM Bugs 15 May 1, 2016 16:12
Problem of simulating of small droplet with radius of 2mm liguifan OpenFOAM Running, Solving & CFD 5 June 3, 2014 02:53
Can anyone help locating an error in my 2-D MATLAB SIMPLE code? Hoggs17 Main CFD Forum 7 July 9, 2013 03:02
[blockMesh] Axisymmetrical mesh Rasmus Gjesing (Gjesing) OpenFOAM Meshing & Mesh Conversion 10 April 2, 2007 14:00


All times are GMT -4. The time now is 16:48.