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

Matlab code for Finite Volume Method in 2D

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

Reply
 
LinkBack Thread Tools Display Modes
Old   November 22, 2014, 02: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!
coagmento is offline   Reply With Quote

Old   December 15, 2014, 05:09
Thumbs up
  #2
New Member
 
Md. Intishar's Avatar
 
Md. Intishar Rahman
Join Date: Nov 2014
Location: Bangladesh
Posts: 6
Rep Power: 2
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   Yesterday, 05:45
Smile Corrected code
  #3
New Member
 
mgedwards
Join Date: Aug 2015
Posts: 1
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
mgedwards is offline   Reply With Quote

Reply

Thread Tools
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 On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
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
alphaEqn.H in twoPhaseEulerFoam cheng1988sjtu OpenFOAM Bugs 14 September 18, 2011 13:46
my new matlab code houkensjtu Main CFD Forum 1 August 26, 2011 00:19
Axisymmetrical mesh Rasmus Gjesing (Gjesing) OpenFOAM Native Meshers: blockMesh 10 April 2, 2007 14:00


All times are GMT -4. The time now is 19:23.