
[Sponsors] 
November 22, 2014, 02:32 
Matlab code for Finite Volume Method in 2D

#1 
New Member
David
Join Date: Nov 2014
Posts: 2
Rep Power: 0 
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 x with initial condition For simplicity and interest, I take , where is the distance function given by so that all the density is concentrated near the point 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((X1/2).^2+(Y1/2).^2); gradphix = (X1/2)./phi; gradphiy = (Y1/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,j1)); %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,j1)); 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(i1,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(i1,j)); end Unew(i,j) = U(i,j)  (delta_t)*(F_rightF_left)(delta_t)*(G_up  ... G_down); end end U = Unew; surf(X,Y,Unew) pause(0.1); end 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 . 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 forloops in my code? My apologies if my question is too elementary for this forum! 

December 15, 2014, 05:09 

#2 
New Member
Md. Intishar Rahman
Join Date: Nov 2014
Location: Bangladesh
Posts: 6
Rep Power: 2 
Hello coagmento,
I have been working on developing a kepsilon 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 

Yesterday, 05:45 
Corrected code

#3 
New Member
mgedwards
Join Date: Aug 2015
Posts: 1
Rep Power: 0 
% 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((X1/2).^2+(Y1/2).^2); gradphix = (X1/2)./phi; gradphiy = (Y1/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,j1)); %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,j1)); 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(i1,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(i1,j)); end Unew(i,j) = U(i,j)  (delta_t)*(F_rightF_left)(delta_t)*(G_up  ... G_down); end end U = Unew; surf(X,Y,Unew) pause(0.1); end 

Thread Tools  
Display Modes  


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 2D 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 