Twelve steps to solve the navier stokes
Twelve steps to solve the navier stokes
Embedded Files
The MATLAB based solvers presented here are based on the works of Prof. Lorena A. Barba, which is originally done in python (Find the work here). A detailed explanation of each step is not included here, instead matlab codes and sample results are provided.
Published on: 25th April, 2023
i. 1-D linear convection
i. 1-D linear convection
Matlab code
length = 2;
nx = 41; %Grid Size in row-direction
deltaX = length/(nx-1);
x = 0:deltaX:length;
nt = 25; %Number of time steps (simulation time)
deltaT = 0.025; %Time Step Size
c = 1; %Wave Speed
u = ones(1,nx);
%Initial Condition
a = floor(0.5/deltaX);
b = floor(1+ (1/deltaX));
u(a:b) = 2;
u_old = u;
%For Animation
p = plot(nan,nan);
p.XData = x;
p.YData = u;
axis manual
xlabel('x')
ylabel('f(x)')
title('Step 1: One Dimensional Linear Convection Equation')
exportgraphics(gca,'1stStep.gif','Append',true)
%Advancement in Time
for t = 1:nt
for i = 2:nx
u(i) = u_old(i) - c*(deltaT/deltaX)*(u_old(i)-u_old(i-1));
end
u_old=u;
p.YData = u;
axis manual
xlabel('x')
ylabel('f(x)')
title('Step 1: 1-D Linear Convection Equation')
exportgraphics(gca,'1stStep.gif','Append',true)
pause(0.10)
end
ii: 1-D non-linear convection
ii: 1-D non-linear convection
Matlab code
length = 2;
nx = 20; %Grid Size in row-direction
deltaX = length/(nx-1);
x = 0:deltaX:length;
nt = 25; %Number of time steps (simulation time)
deltaT = 0.025; %Time Step Size
u = ones(1,nx);
%Initial Condition
a = floor(0.5/deltaX);
b = floor(1+ (1/deltaX));
u(a:b) = 2;
u_old = u;
%For Animation
p = plot(nan,nan);
p.XData = x;
p.YData = u;
axis manual
xlabel('x')
ylabel('f(x)')
title('Step 2: 1-D Non-linear Convection Equation')
exportgraphics(gcf,'2ndStep.gif','Append',true)
%Advancement in Time
for t = 1:nt
for i = 2:nx
u(i) = u_old(i) - u_old(i)*(deltaT/deltaX)*(u_old(i)-u_old(i-1));
end
u_old=u;
p.YData = u;
axis manual
xlabel('x')
ylabel('f(x)')
title('Step 2: 1-D Non-linear Convection Equation')
exportgraphics(gcf,'2ndStep.gif','Append',true)
pause(0.10)
end
iii: 1-D diffusion equation
iii: 1-D diffusion equation
Matlab code
length = 2;
nx = 40; %Grid Size in row-direction
deltaX = length/(nx-1);
x = 0:deltaX:length;
nt = 25; %Number of time steps (simulation time)
sigma = 0.2; % Parameter
nu = 0.3; % Viscosity
u = ones(1,nx);
deltaT = (sigma*deltaX*deltaX/nu); % Time Step Size
u = ones(1,nx);
%Initial Condition
a= floor(0.5/deltaX);
b= floor(1+ (1/deltaX));
u(a:b) = 2;
u_old = u;
%For Animation
p = plot(nan,nan);
p.XData = x;
p.YData = u;
axis manual
xlabel('x')
ylabel('f(x)')
title('Step 3: 1-D Diffusion Equation')
exportgraphics(gcf,'3rdStep.gif','Append',true)
%Advancement in Time
for t = 1:nt
for i = 2:nx-1
u(i) = u_old(i) + nu*(deltaT/(deltaX^2))*(u_old(i+1)-2*u_old(i)+u_old(i-1));
end
u_old=u;
p.YData = u;
axis manual
xlabel('x')
ylabel('f(x)')
title('Step 3: 1-D Diffusion Equation')
exportgraphics(gcf,'3rdStep.gif','Append',true)
pause(0.10)
end
iv: 1-D burger's equation
iv: 1-D burger's equation
Matlab code
length = 2*pi;
nx = 40; %Grid Size in row-direction
deltaX = length/(nx-1);
nt = 20; %Number of time steps (simulation time)
sigma = 0.2; % Parameter
nu = 0.07; % Viscosity
%deltaT = (sigma*deltaX*deltaX/nu); % Time Step Size
deltaT = deltaX*nu; % Time Step Size
u = ones(1,nx);
u_old = ones(1,nx);
u_analytical = ones(1,nx);
xVal = 0:deltaX:length;
%Initial Condition
%for i>1
xVal_ = xVal;
nu_ = nu;
syms xVal_ nu_ t_
phi = exp(-(xVal_ - 4 * t_)^2 / (4 * nu_ * (t_ + 1))) + exp(- (xVal_ - 4 * t_ - 2 * sym(pi))^2 / (4 * nu_ * (t_ + 1)));
dphix = diff(phi,xVal_);
fn = -2*nu_*dphix/phi +4 ;
fn1 = subs(fn,t_,0);
fn2 = subs(fn1,nu_,0.07);
xrange = linspace(0,2*pi,nx);
u = eval(subs(fn2,xVal_,xVal));
u_old = u;
%For Animation
p = plot(nan,nan);
p.XData = xVal;
p.YData = u;
axis manual
xlabel('x')
ylabel('f(x)')
title('Step 4: 1-D Burgers Equation')
exportgraphics(gcf,'4thStep.gif','Append',true)
for t = 1:nt
for i = 2:nx-1
u(i) = u_old(i) - u_old(i)*(deltaT/deltaX)*(u_old(i)-u_old(i-1)) + nu*(deltaT/(deltaX*deltaX))*(u_old(i+1)-2*u_old(i)+u_old(i-1));
end
%Periodic Boundary Condition
u(nx) = u_old(nx) - u_old(nx)*(deltaT/deltaX)*(u_old(nx)-u_old(nx-1)) + nu*(deltaT/(deltaX*deltaX))*(u_old(1)-2*u_old(nx)+u_old(nx-1));
u(1) = u_old(1) - u_old(1)*(deltaT/deltaX)*(u_old(1)-u_old(nx)) + nu*(deltaT/(deltaX*deltaX))*(u_old(2)-2*u_old(1)+u_old(nx));
u_old = u;
p.YData = u;
axis manual
xlabel('x')
ylabel('f(x)')
title('Step 4: 1-D Burgers Equation')
exportgraphics(gcf,'4thStep.gif','Append',true)
pause(0.10)
end
fn1 = subs(fn,t_,(nt*deltaT));
fn2 = subs(fn1,nu_,0.07);
u_analytical = eval(subs(fn2,xVal_,xVal));
yyaxis left
plot(xVal,u,'r');
hold on
plot(xVal,u_analytical,'b--');25
v: 2-D linear convection
v: 2-D linear convection
Matlab code
length = 2; %Domain Length
height = 2; %Domain Height
nx = 81; %Grid Size in x- direction
ny = 81; %Grid Size in y- direction
deltaX = length/(nx-1);
deltaY = height/(ny-1);
nt = 100; %Number of time steps (simulation time)
sigma = 0.25;
deltaT = sigma*min(deltaX,deltaY);
x = (0:deltaX:length);
y = (0:deltaY:height);
c = 1; %Wave Speed
%Initialize u
u = ones(ny, nx);
%Intitial Conditions for u
a = floor(0.5/deltaX);
b = floor(1+ (1/deltaX));
p = floor(0.5/deltaY);
q = floor(1+ (1/deltaY));
u(p:q,a:b) = 2;
%Boundary Conditions
u(1:end,end) = 1;
u(1:end,end) = 1;
u(1,2:end-1) = 1;
u(end,2:end-1)= 1;
u_old = u;
[X,Y] =meshgrid(x, y) ;
colormap(jet)
subplot(1,2,1)
surf(X,Y,u)
xlabel('x')
ylabel('y')
zlabel('f(x,y)')
title('Step 5: 2-D Linear Convection Equation')
subplot(1,2,2)
contour(X,Y,u)
axis tight manual
xlabel('x')
ylabel('y')
exportgraphics(gcf,'5thStep.gif','Append',true)
%set(gca,'nextplot','replacechildren');
%v = VideoWriter('step5.avi');
%open(v);
%Time Marching
for t = 1:nt
for j = 2:ny-1
for i = 2:nx-1
u(j,i) = u_old(j,i) - c*(deltaT/deltaX)*(u_old(j,i)-u_old(j,i-1))- c*(deltaT/deltaY)*(u_old(j,i)-u_old(j-1,i));
end
end
u_old=u;
colormap(jet)
subplot(1,2,1)
surf(X,Y,u)
xlabel('x')
ylabel('y')
zlabel('f(x,y)')
title('Step 5: 2-D Linear Convection Equation')
subplot(1,2,2)
contour(X,Y,u)
xlabel('x')
ylabel('y')
exportgraphics(gcf,'5thStep.gif','Append',true)
%frame = getframe(gcf);
%writeVideo(v,frame);
pause(0.10)
end
vi: 2-D non-linear convection
vi: 2-D non-linear convection
Matlab code
length = 2; %Domain Length
height = 2; %Domain Height
nx = 81; %Grid Size in x- direction
ny = 81; %Grid Size in y- direction
deltaX = length/(nx-1);
deltaY = height/(ny-1);
nt = 100; %Number of time steps (simulation time)
sigma = 0.25;
deltaT = sigma*min(deltaX,deltaY);
x = (0:deltaX:length);
y = (0:deltaY:height);
c = 1; %Wave Speed
%FunctionInitialization
u = ones(ny,nx);
v = ones(ny,nx);
%Intitial Conditions for u
a = floor(0.5/deltaX);
b = floor(1+ (1/deltaX));
p = floor(0.5/deltaY);
q = floor(1+ (1/deltaY));
u(a:b, p:q) = 2.0;
v(a:b, p:q) = 2.0;
%Boundary Conditions
u(1:end,1) = 1;
u(1:end,end) = 1;
v(1:end,1) = 1;
v(1:end,end) = 1;
%top&bottom
u(1,2:end-1) = 1;
u(end,2:end-1)= 1;
v(1,2:end-1) = 1;
v(end,2:end-1)= 1;
%%%%%%%%%%%%%%%%%%%
u_old = u;
v_old = v;
[X,Y] =meshgrid(x, y) ;
colormap(jet)
subplot(1,2,1)
surf(X,Y,u)
xlabel('x')
ylabel('y')
zlabel('f(x,y)')
title('Step 6: 2-D Non-Linear Convection Equation')
subplot(1,2,2)
contour(X,Y,u)
xlabel('x')
ylabel('y')
axis tight manual
exportgraphics(gcf,'6thStep.gif','Append',true)
%set(gca,'nextplot','replacechildren');
%v = VideoWriter('step5.avi');
%open(v);
%Time Marching
for t = 1:nt
for j = 2:ny-1
for i = 2:nx-1
u(j,i) = u_old(j,i) - u_old(j,i)*(deltaT/deltaX)*(u_old(j,i)-u_old(j,i-1))- v_old(j,i)*(deltaT/deltaY)*(u_old(j,i)-u_old(j-1,i));
v(j,i) = v_old(j,i) - u_old(j,i)*(deltaT/deltaX)*(v_old(j,i)-v_old(j,i-1))- v_old(j,i)*(deltaT/deltaY)*(v_old(j,i)-v_old(j-1,i));
end
end
u_old=u;
v_old=v;
colormap(jet)
subplot(1,2,1)
surf(X,Y,u)
xlabel('x')
ylabel('y')
zlabel('f(x,y)')
title('Step 6: 2-D Non-Linear Convection Equation')
subplot(1,2,2)
contour(X,Y,u)
xlabel('x')
ylabel('y')
axis tight manual
exportgraphics(gcf,'6thStep.gif','Append',true)
%frame = getframe(gcf);
%writeVideo(v,frame);
pause(0.10)
end
vii: 2-D diffusion
vii: 2-D diffusion
Matlab code
length = 2; %Domain Length
height = 2; %Domain Height
nx = 31; %Grid Size in x- direction
ny = 31; %Grid Size in y- direction
deltaX = length/(nx-1);
deltaY = height/(ny-1);
nt = 30; %Number of time steps (simulation time)
sigma = 0.25;
nu = 0.05;
deltaT = (sigma*deltaX*deltaY)/nu;
x = (0:deltaX:length);
y = (0:deltaY:height);
%FunctionInitialization
u = ones(ny,nx);
v = ones(ny,nx);
%Intitial Conditions for u
a = floor(0.5/deltaX);
b = floor(1+ (1/deltaX));
p = floor(0.5/deltaY);
q = floor(1+ (1/deltaY));
u(a:b, p:q) = 2.0;
v(a:b, p:q) = 2.0;
%Boundary Conditions
u(1:end,1) = 1;
u(1:end,end) = 1;
v(1:end,1) = 1;
v(1:end,end) = 1;
%top&bottom
u(1,2:end-1) = 1;
u(end,2:end-1)= 1;
v(1,2:end-1) = 1;
v(end,2:end-1)= 1;
%%%%%%%%%%%%%%%%%%%
u_old = u;
v_old = v;
%%%%%%%%%%%%%%%%%%%%%%%%%%
[X,Y] =meshgrid(x, y) ;
colormap(jet)
subplot(1,2,1)
surf(X,Y,u)
xlabel('x')
ylabel('y')
zlabel('f(x,y)')
title('Step 7: 2-D Diffusion Equation')
subplot(1,2,2)
contour(X,Y,u)
xlabel('x')
ylabel('y')
axis tight manual
exportgraphics(gcf,'7thStep.gif','Append',true)
%set(gca,'nextplot','replacechildren');
%v = VideoWriter('step5.avi');
%open(v);
%%%%%%%%%%%%%%%%%%%%%%%%%%
%Time Marching
for t = 1:nt
for j = 2:ny-1
for i = 2:nx-1
u(j,i) = u_old(j,i) + nu*(deltaT/(deltaX*deltaX))*(u_old(j,i+1)-2*u_old(j,i)+u_old(j,i-1)) + nu*(deltaT/(deltaY*deltaY))*(u_old(j+1,i)-2*u_old(j,i)+u_old(j-1,i));
v(j,i) = v_old(j,i) + nu*(deltaT/(deltaX*deltaX))*(v_old(j,i+1)-2*v_old(j,i)+v_old(j,i-1)) + nu*(deltaT/(deltaY*deltaY))*(v_old(j+1,i)-2*v_old(j,i)+v_old(j-1,i));
end
end
u_old=u;
v_old=v;
colormap(jet)
subplot(1,2,1)
surf(X,Y,u)
xlabel('x')
ylabel('y')
zlabel('f(x,y)')
title('Step 7: 2-D Diffusion Equation')
subplot(1,2,2)
contour(X,Y,u)
xlabel('x')
ylabel('y')
axis tight manual
exportgraphics(gcf,'7thStep.gif','Append',true)
pause(0.10)
end
viii: 2-D diffusion
viii: 2-D diffusion
Matlab code
length = 2; %Domain Length
height = 2; %Domain Height
nx = 41; %Grid Size in x- direction
ny = 41; %Grid Size in y- direction
deltaX = length/(nx-1);
deltaY = height/(ny-1);
nt = 120; %Number of time steps (simulation time)
sigma = 0.0009;
nu = 0.01;
c = 1.0; %Wave Velocity
deltaT = (sigma*deltaX*deltaY)/nu;
x = (0:deltaX:length);
y = (0:deltaY:height);
%FunctionInitialization
u = ones(ny,nx);
v = ones(ny,nx);
%Intitial Conditions for u
a = floor(0.5/deltaX);
b = floor(1+ (1/deltaX));
p = floor(0.5/deltaY);
q = floor(1+ (1/deltaY));
u(a:b, p:q) = 2.0;
v(a:b, p:q) = 2.0;
%Boundary Conditions
u(1:end,1) = 1;
u(1:end,end) = 1;
v(1:end,1) = 1;
v(1:end,end) = 1;
%top&bottom
u(1,2:end-1) = 1;
u(end,2:end-1)= 1;
v(1,2:end-1) = 1;
v(end,2:end-1)= 1;
%%%%%%%%%%%%%%%%%%%
u_old = u;
v_old = v;
%%%%%%%%%%%%%%%%%%%%%%%%%%
%___Create GIF___
[X,Y] =meshgrid(x, y) ;
colormap(jet)
subplot(1,2,1)
surf(X,Y,u)
xlabel('x')
ylabel('y')
zlabel('f(x,y)')
title('Step 8: 2-D Burgers Equation')
subplot(1,2,2)
contour(X,Y,u)
xlabel('x')
ylabel('y')
axis tight manual
exportgraphics(gcf,'8thStep.gif','Append',true)
%set(gca,'nextplot','replacechildren');
%v = VideoWriter('step5.avi');
%open(v);
%%%%%%%%%%%%%%%%%%%%%%%%%%
%Time Marching
for t = 1:nt
for j = 2:ny-1
for i = 2:nx-1
u(j,i) = u_old(j,i) -(deltaT/deltaX)*u_old(j,i)*(u_old(j,i)-u_old(j,i-1)) - (deltaT/deltaY)*v_old(j,i)*(u_old(j,i)-u_old(j-1,i)) + nu*(deltaT/(deltaX*deltaX))*(u_old(j,i+1)-2*u_old(j,i)+u_old(j,i-1)) + nu*(deltaT/(deltaY*deltaY))*(u_old(j+1,i)-2*u_old(j,i)+u_old(j-1,i));
v(j,i) = v_old(j,i) -(deltaT/deltaX)*u_old(j,i)*(v_old(j,i)-v_old(j,i-1)) - (deltaT/deltaY)*v_old(j,i)*(v_old(j,i)-v_old(j-1,i)) + nu*(deltaT/(deltaX*deltaX))*(v_old(j,i+1)-2*v_old(j,i)+v_old(j,i-1)) + nu*(deltaT/(deltaY*deltaY))*(v_old(j+1,i)-2*v_old(j,i)+v_old(j-1,i));
end
end
u_old=u;
v_old=v;
%%%%%%%%%%%%%%%%%%%%%%%%%%
%___Create GIF___
colormap(jet)
subplot(1,2,1)
surf(X,Y,u)
xlabel('x')
ylabel('y')
zlabel('f(x,y)')
title('Step 8: 2-D Burgers Equation')
subplot(1,2,2)
contour(X,Y,u)
xlabel('x')
ylabel('y')
axis tight manual
exportgraphics(gcf,'8thStep.gif','Append',true)
%set(gca,'nextplot','replacechildren');
%v = VideoWriter('step5.avi');
%open(v);
pause(0.10)
%%%%%%%%%%%%%%%%%%%%%%%%%%
end
ix: 2-D laplace equation
ix: 2-D laplace equation
Matlab code
length = 2; %Domain Length
height = 1; %Domain Heightclc
nx = 31; %Grid Size in x- direction
ny = 31; %Grid Size in y- direction
deltaX = length/(nx-1);
deltaY = height/(ny-1);
x = (0:deltaX:length);
y = (0:deltaY:height);
tolerance = 1e-4;
%Initial Condition
p = zeros(ny, nx);
%Boundary Conditions
% @ X = 0
p(1,:) = 0;
% @ X = Length
p(1:ny,nx) = y(1:ny);
% @ y=0 (for Interior points in x-direction)
p(1,2:nx-1) = p(2,2:nx-1);
% @ y=height (for Interior points in x-direction)
p(ny,2:nx-1) = p(ny-1,2:nx-1);
%%%%%%%%%%%%%%%%%%%%%%%%%%
%___Create GIF___
[X,Y] =meshgrid(x, y) ;
colormap(jet)
%subplot(1,2,1)
surf(X,Y,p)
%subplot(1,2,2)
%contour(X,Y,p)
xlabel('x')
ylabel('y')
zlabel('f(x,y)')
title('Step 9: 2-D Laplace Equation')
axis tight manual
exportgraphics(gcf,'9thStep.gif','Append',true)
%set(gca,'nextplot','replacechildren');
%v = VideoWriter('step5.avi');
%open(v);
%%%%%%%%%%%%%%%%%%%%%%%%%%
count = 0;
p_norm = 1.0;
while p_norm>=tolerance
p_old = p;
for i = 2:nx-1
for j = 2:ny-1
p(j,i) = ( (deltaY*deltaY)*(p(j,i+1)+p(j,i-1)) + (deltaX*deltaX)*(p(j+1,i)+p(j-1,i)) ) / (2.0*(deltaX*deltaX + deltaY*deltaY));
end
end
%Boundary Conditions
% @ X = 0
p(1,:) = 0;
% @ X = Length
p(1:ny,nx) = y(1:ny);
% @ y=0 (for Interior points in x-direction)
p(1,2:nx-1) = p(2,2:nx-1);
% @ y=height (for Interior points in x-direction)
p(ny,2:nx-1) = p(ny-1,2:nx-1);
p_norm = abs(norm(p)-norm(p_old))/norm(abs(p_old));
count = count + 1;
if(mod(count,20)==0)
%%%%%%%%%%%%%%%%%%%%%%%%%%
%___Create GIF___
[X,Y] =meshgrid(x, y) ;
colormap(jet)
%subplot(1,2,1)
surf(X,Y,p)
%subplot(1,2,2)
%contour(X,Y,p)
xlabel('x')
ylabel('y')
zlabel('f(x,y)')
title('Step 9: 2-D Laplace Equation')
axis tight manual
exportgraphics(gcf,'9thStep.gif','Append',true)
%set(gca,'nextplot','replacechildren');
%v = VideoWriter('step5.avi');
%open(v);
pause(0.01)
%%%%%%%%%%%%%%%%%%%%%%%%%%
end
end
fprintf("Laplace Equation Converged at %d steps\n",count);
x: 2-D poisson equation
x: 2-D poisson equation
Matlab code
length = 2; %Domain Length
height = 1; %Domain Height
nx = 30; %Grid Size in x- direction
ny = 30; %Grid Size in y- direction
nt = 100;
deltaX = length/(nx-1);
deltaY = height/(ny-1);
x = (0:deltaX:length);
y = (0:deltaY:height);
tolerance = 1e-4;
%Initial Condition
p = zeros(ny,nx);
b = zeros(ny,nx);
%Boundary Conditions
% @ X = 0
p(1,:) = 0;
% @ X = Length
p(1:ny,nx) = 0;
% @ y=0 (for Interior points in x-direction)
p(1,2:nx-1) = 0;
% @ y=height (for Interior points in x-direction)
p(ny,2:nx-1) = 0;
b(ceil(ny/4),ceil(nx/4)) = 100;
b(ceil((3*ny)/4),ceil((3*nx)/4)) = -100;
%%%%%%%%%%%%%%%%%%%%%%%%%%
%___Create GIF___
[X,Y] =meshgrid(x, y) ;
colormap(jet)
%subplot(1,2,1)
surf(X,Y,p)
%subplot(1,2,2)
%contour(X,Y,p)
xlabel('x')
ylabel('y')
zlabel('f(x,y)')
title('Step 10: 2-D Poisson Equation')
axis tight manual
exportgraphics(gcf,'10thStep.gif','Append',true)
%set(gca,'nextplot','replacechildren');
%v = VideoWriter('step5.avi');
%open(v);
%%%%%%%%%%%%%%%%%%%%%%%%%%
count = 0;
p_norm = 1.0;
%while p_norm>tolerance
while count<=nt
p_old = p;
for i = 2:nx-1
for j = 2:ny-1
rhs = b(j,i)*deltaX*deltaX*deltaY*deltaY;
p(j,i) = ( (deltaY*deltaY)*(p(j,i+1)+p(j,i-1)) + (deltaX*deltaX)*(p(j+1,i)+p(j-1,i)) - rhs) / (2.0*(deltaX*deltaX + deltaY*deltaY));
end
end
%Boundary Conditions
% @ X = 0
p(1,:) = 0;
% @ X = Length
p(1:ny,nx) = 0;
% @ y=0 (for Interior points in x-direction)
p(1,2:nx-1) = 0;
% @ y=height (for Interior points in x-direction)
p(ny,2:nx-1) = 0;
b(ceil(ny/4),ceil(nx/4)) = 100;
b(ceil((3*ny)/4),ceil((3*nx)/4)) = -100;
p_norm = abs(norm(p)-norm(p_old))/norm(abs(p_old));
if(mod(count,10)==0)
%%%%%%%%%%%%%%%%%%%%%%%%%%
%___Create GIF___
[X,Y] =meshgrid(x, y) ;
colormap(jet)
%subplot(1,2,1)
surf(X,Y,p)
%subplot(1,2,2)
%contour(X,Y,p)
xlabel('x')
ylabel('y')
zlabel('f(x,y)')
title('Step 10: 2-D Poisson Equation')
axis tight manual
exportgraphics(gcf,'10thStep.gif','Append',true)
%set(gca,'nextplot','replacechildren');
%v = VideoWriter('step5.avi');
%open(v);
pause(0.01)
%%%%%%%%%%%%%%%%%%%%%%%%%%
end
count = count + 1;
end
p_norm
fprintf("Laplace Equation Converged at %d steps\n",count);
xi: 2-D cavity problem
xi: 2-D cavity problem
Matlab code
length = 2; %Domain Length
height = 2; %Domain Height
global nx ny nt nu deltaX deltaY deltaT x y rho tolerance ppe_count;
ppe_count = 1000;
nx = 50; %Grid Size in x- direction
ny = 50; %Grid Size in y- direction
nt = 1000;
nu = 0.1;
rho = 1.0;
deltaX = length/(nx-1);
deltaY = height/(ny-1);
deltaT = 0.001;
x = (0:deltaX:length);
y = (0:deltaY:height);
tolerance = 1e-4;
time= (0:deltaT:(nt*deltaT));
p_norm = 1;
count = 0;
%Initial Condition
p = zeros(ny,nx);
b = zeros(ny,nx);
u = zeros(ny,nx);
v = zeros(ny,nx);
V_magn = zeros(ny,nx);
%Boundary Conditions
u(:,1) = 0;
u(:,nx) = 0;
v(:,1) = 0;
v(:,nx) = 0;
u(1, 2:nx-1) = 0;
u(ny, 2:nx-1) = 1;
v(1, 2:nx-1) = 0;
v(ny, 2:nx-1) = 0;
p(:,1) = p(:,2);
p(:,nx) = p(:,nx-1);
p(1,:) = p(2,:);
p(ny,:) = 0;
V_magn(:,:) = sqrt(u(:,:).*u(:,:) + v(:,:).*v(:,:));
%%%%%%%%%%%%%%%%%%%%%%%%%%
%___Create GIF___
[X,Y] =meshgrid(x, y) ;
colormap(jet)
%subplot(1,2,1)
%surf(X,Y,V_magn)
%subplot(1,2,2)
contourf(X,Y,V_magn)
xlabel('x')
ylabel('y')
legend('Velocity Magnitude')
title('Step 11: 2-D Cavity Flow')
axis tight manual
exportgraphics(gcf,'11thStep.gif','Append',true)
%set(gca,'nextplot','replacechildren');
%v = VideoWriter('step5.avi');
%open(v);
%%%%%%%%%%%%%%%%%%%%%%%%%%
%Time Marching
for t = 1:nt
u_old = u;
v_old = v;
for j = 2:ny-1
for i = 2:nx-1
convectionTermsX = u_old(j,i)*((u_old(j,i)-u_old(j,i-1))/deltaX) + v_old(j,i)*((u_old(j,i)-u_old(j-1,i))/deltaY);
convectionTermsY = u_old(j,i)*((v_old(j,i)-v_old(j,i-1))/deltaX) + v_old(j,i)*((v_old(j,i)-v_old(j-1,i))/deltaY);
DiffusionX = nu*(u_old(j,i+1)-2*u_old(j,i)+u_old(j,i-1))/(deltaX*deltaX) + nu*(u_old(j+1,i)-2*u_old(j,i)+u_old(j-1,i))/(deltaY*deltaY);
DiffusionY = nu*(v_old(j,i+1)-2*v_old(j,i)+v_old(j,i-1))/(deltaX*deltaX) + nu*(v_old(j+1,i)-2*v_old(j,i)+v_old(j-1,i))/(deltaY*deltaY);
PressureX = (p(j,i+1)-p(j,i-1))/(rho*2*deltaX);
PressureY = (p(j+1,i)-p(j-1,i))/(rho*2*deltaY);
u(j,i) = u_old(j,i) + deltaT*(DiffusionX - PressureX- convectionTermsX);
v(j,i) = v_old(j,i) + deltaT*(DiffusionY - PressureY- convectionTermsY);
end
end
%Boundary Conditions
u(:,1) = 0;
u(:,nx) = 0;
v(:,1) = 0;
v(:,nx) = 0;
u(1, 2:nx-1) = 0;
u(ny, 2:nx-1) = 1;
v(1, 2:nx-1) = 0;
v(ny, 2:nx-1) = 0;
p(:,1) = p(:,2);
p(:,nx) = p(:,nx-1);
p(1,:) = p(2,:);
p(ny,:) = 0;
[p, count, p_norm] = ppe_cavity(p, u, v);
if(mod(t,25) == 0 || t ==1 )
dummy1 = mean(u, 'all');
dummy2 = mean(v, 'all');
V_avg = sqrt(dummy1*dummy1 + dummy2*dummy2);
% fprintf("Avg Velocity is = %f\n", V_avg);
plot(time(t),V_avg,'b--o')
% fprintf("The Press-Poisson Eqn converged at %d and with tolerance = %f\n", count, p_norm);
V_magn(:,:) = sqrt(u(:,:).*u(:,:) + v(:,:).*v(:,:));
%%%%%%%%%%%%%%%%%%%%%%%%%%
%___Create GIF___
colormap(jet)
%subplot(1,2,1)
%surf(X,Y,V_magn)
%subplot(1,2,2)
contourf(X,Y,V_magn)
xlabel('x')
ylabel('y')
legend('Velocity Magnitude')
title('Step 11: 2-D Cavity Flow')
axis tight manual
exportgraphics(gcf,'11thStep.gif','Append',true)
%set(gca,'nextplot','replacechildren');
%v = VideoWriter('step5.avi');
%open(v);
pause(0.10);
%%%%%%%%%%%%%%%%%%%%%%%%%%
% fprintf("Poisson Equation converged at %d\n",count);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[p_new, count, p_norm] = ppe_cavity(p, u, v)
count = 0;
p_norm = 1.0;
p_new = p;
global nx ny deltaX deltaY deltaT rho ppe_count;
while count < ppe_count
% p_norm < tolerance
p_old = p_new;
for i = 2:nx-1
for j = 2:ny-1
dudx = (u(j,i+1)-u(j,i-1))/(2*deltaX);
dvdx = (v(j,i+1)-v(j,i-1))/(2*deltaX);
dudy = (u(j+1,i)-u(j-1,i))/(2*deltaY);
dvdy = (v(j+1,i)-v(j-1,i))/(2*deltaY);
b = ((1/deltaT)*(dudx + dvdy))-(dudx*dudx + 2*dudy*dvdx + dvdy*dvdy);
rhs = b*rho*deltaX*deltaX*deltaY*deltaY;
p_new(j,i) = ( (deltaY*deltaY)*(p(j,i+1)+p(j,i-1)) + (deltaX*deltaX)*(p(j+1,i)+p(j-1,i)) - rhs) / (2.0*(deltaX*deltaX + deltaY*deltaY));
end
end
%Boundary Conditions
p_new(:,1) = p_new(:,2);
p_new(:,nx) = p_new(:,nx-1);
p_new(1,:) = p_new(2,:);
p_new(ny,:) = 0;
p_norm = abs(norm(p_new)-norm(p_old))/norm(abs(p_old));
count = count + 1;
end
end
xii: 2-D channel problem
xii: 2-D channel problem
Matlab code
length = 2; %Domain Length
height = 2; %Domain Height
global nx ny nt nu deltaX deltaY deltaT x y rho tolerance ppe_count;
ppe_count = 50;
nx = 41; %Grid Size in x- direction
ny = 41; %Grid Size in y- direction
nt = 300;
nu = 0.1;
F = 1.0;
rho = 1.0;
deltaX = length/(nx-1);
deltaY = height/(ny-1);
deltaT = 0.01;
x = (0:deltaX:length);
y = (0:deltaY:height);
tolerance = 1e-4;
time= (0:deltaT:(nt*deltaT));
p_norm = 1;
count = 0;
%Initial Condition
p = zeros(ny,nx);
b = zeros(ny,nx);
u = zeros(ny,nx);
v = zeros(ny,nx);
V_magn = zeros(ny,nx);
%Boundary Conditions
u(1,:) = 0;
u(ny,:) = 0;
v(1,:) = 0;
v(ny,:) = 0;
p(1,:) = p(2,:);
p(ny,:) = p(ny-1,:);
V_magn(:,:) = sqrt(u(:,:).*u(:,:) + v(:,:).*v(:,:));
%%%%%%%%%%%%%%%%%%%%%%%%%%
%___Create GIF___
[X,Y] =meshgrid(x, y) ;
%For_GIF
%colormap(jet)
%contourf(X,Y,V_magn)
%axis tight manual
%exportgraphics(gcf,'12thStep.gif','Append',true)
%For_Video
%set(gca,'nextplot','replacechildren');
%v = VideoWriter('step5.avi');
%open(v);
%%%%%%%%%%%%%%%%%%%%%%%%%%
%Time Marching
for t = 1:nt
u_old = u;
v_old = v;
for j = 2:ny-1
for i = 2:nx-1
convectionTermsX = u_old(j,i)*((u_old(j,i)-u_old(j,i-1))/deltaX) + v_old(j,i)*((u_old(j,i)-u_old(j-1,i))/deltaY);
convectionTermsY = u_old(j,i)*((v_old(j,i)-v_old(j,i-1))/deltaX) + v_old(j,i)*((v_old(j,i)-v_old(j-1,i))/deltaY);
DiffusionX = nu*(u_old(j,i+1)-2*u_old(j,i)+u_old(j,i-1))/(deltaX*deltaX) + nu*(u_old(j+1,i)-2*u_old(j,i)+u_old(j-1,i))/(deltaY*deltaY);
DiffusionY = nu*(v_old(j,i+1)-2*v_old(j,i)+v_old(j,i-1))/(deltaX*deltaX) + nu*(v_old(j+1,i)-2*v_old(j,i)+v_old(j-1,i))/(deltaY*deltaY);
PressureX = (p(j,i+1)-p(j,i-1))/(rho*2*deltaX);
PressureY = (p(j+1,i)-p(j-1,i))/(rho*2*deltaY);
u(j,i) = u_old(j,i) + deltaT*(DiffusionX - PressureX- convectionTermsX + F);
v(j,i) = v_old(j,i) + deltaT*(DiffusionY - PressureY- convectionTermsY);
end
end
%Boundary Conditions
u(1,:) = 0;
u(ny,:) = 0;
v(1,:) = 0;
v(ny,:) = 0;
p(1,:) = p(2,:);
p(ny,:) = p(ny-1,:);
for j = 2:ny-1
% @ i = 1
i = 1;
convectionTermsX = u_old(j,i)*((u_old(j,i)-u_old(j,nx))/deltaX) + v_old(j,i)*((u_old(j,i)-u_old(j-1,i))/deltaY);
convectionTermsY = u_old(j,i)*((v_old(j,i)-v_old(j,nx))/deltaX) + v_old(j,i)*((v_old(j,i)-v_old(j-1,i))/deltaY);
DiffusionX = nu*(u_old(j,i+1)-2*u_old(j,i)+u_old(j,nx))/(deltaX*deltaX) + nu*(u_old(j+1,i)-2*u_old(j,i)+u_old(j-1,i))/(deltaY*deltaY);
DiffusionY = nu*(v_old(j,i+1)-2*v_old(j,i)+v_old(j,nx))/(deltaX*deltaX) + nu*(v_old(j+1,i)-2*v_old(j,i)+v_old(j-1,i))/(deltaY*deltaY);
PressureX = (p(j,i+1)-p(j,nx))/(rho*2*deltaX);
PressureY = (p(j+1,i)-p(j-1,i))/(rho*2*deltaY);
u(j,i) = u_old(j,i) + deltaT*(DiffusionX - PressureX- convectionTermsX + F);
v(j,i) = v_old(j,i) + deltaT*(DiffusionY - PressureY- convectionTermsY);
% @ i = nx
i = nx;
convectionTermsX = u_old(j,i)*((u_old(j,i)-u_old(j,i-1))/deltaX) + v_old(j,i)*((u_old(j,i)-u_old(j-1,i))/deltaY);
convectionTermsY = u_old(j,i)*((v_old(j,i)-v_old(j,i-1))/deltaX) + v_old(j,i)*((v_old(j,i)-v_old(j-1,i))/deltaY);
DiffusionX = nu*(u_old(j,1)-2*u_old(j,i)+u_old(j,i-1))/(deltaX*deltaX) + nu*(u_old(j+1,i)-2*u_old(j,i)+u_old(j-1,i))/(deltaY*deltaY);
DiffusionY = nu*(v_old(j,1)-2*v_old(j,i)+v_old(j,i-1))/(deltaX*deltaX) + nu*(v_old(j+1,i)-2*v_old(j,i)+v_old(j-1,i))/(deltaY*deltaY);
PressureX = (p(j,1)-p(j,i-1))/(rho*2*deltaX);
PressureY = (p(j+1,i)-p(j-1,i))/(rho*2*deltaY);
u(j,i) = u_old(j,i) + deltaT*(DiffusionX - PressureX- convectionTermsX + F);
v(j,i) = v_old(j,i) + deltaT*(DiffusionY - PressureY- convectionTermsY);
end
[p, count, p_norm] = ppe_channel(p, u, v);
if(mod(t,25) == 0 || t ==1 )
dummy1 = mean(u, 'all');
dummy2 = mean(v, 'all');
V_avg = sqrt(dummy1*dummy1 + dummy2*dummy2);
% fprintf("Avg Velocity is = %f\n", V_avg);
%plot(time(t),V_avg,'b--o')
%fprintf("The Press-Poisson Eqn converged at %d and with tolerance = %f\n", count, p_norm);
V_magn(:,:) = sqrt(u(:,:).*u(:,:) + v(:,:).*v(:,:));
%%%%%%%%%%%%%%%%%%%%%%%%%%
%___Create GIF___
colormap(jet)
%subplot(1,2,1)
%surf(X,Y,V_magn)
%subplot(1,2,2)
contourf(X,Y,V_magn)
xlabel('x')
ylabel('y')
%legend('Velocity Magnitude')
title('Step 12: 2-D Channel Flow')
axis tight manual
exportgraphics(gcf,'12thStep.gif','Append',true)
%set(gca,'nextplot','replacechildren');
%v = VideoWriter('step5.avi');
%open(v);
pause(1.5)
%%%%%%%%%%%%%%%%%%%%%%%%%%
% fprintf("Poisson Equation converged at %d\n",count);
end
end
quiver(x,y,u,v,2)
xlabel('x')
ylabel('y')
title('Step 12: 2-D Channel Flow')
axis equal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[p_new, count, p_norm] = ppe_channel(p, u, v)
count = 0;
p_norm = 1.0;
p_new = p;
global nx ny deltaX deltaY deltaT rho ppe_count;
while count < ppe_count
% p_norm < tolerance
p_old = p_new;
for i = 2:nx-1
for j = 2:ny-1
dudx = (u(j,i+1)-u(j,i-1))/(2*deltaX);
dvdx = (v(j,i+1)-v(j,i-1))/(2*deltaX);
dudy = (u(j+1,i)-u(j-1,i))/(2*deltaY);
dvdy = (v(j+1,i)-v(j-1,i))/(2*deltaY);
b = ((1/deltaT)*(dudx + dvdy))-(dudx*dudx + 2*dudy*dvdx + dvdy*dvdy);
rhs = b*rho*deltaX*deltaX*deltaY*deltaY;
p_new(j,i) = ( (deltaY*deltaY)*(p(j,i+1)+p(j,i-1)) + (deltaX*deltaX)*(p(j+1,i)+p(j-1,i)) - rhs) / (2.0*(deltaX*deltaX + deltaY*deltaY));
end
end
%Boundary Conditions
p_new(1,:) = p_new(2,:);
p_new(ny,:) = p_new(ny-1,:);
for j = 2:ny-1
i = nx;
dudx = (u(j,1)-u(j,i-1))/(2*deltaX);
dvdx = (v(j,1)-v(j,i-1))/(2*deltaX);
dudy = (u(j+1,i)-u(j-1,i))/(2*deltaY);
dvdy = (v(j+1,i)-v(j-1,i))/(2*deltaY);
b = ((1/deltaT)*(dudx + dvdy))-(dudx*dudx + 2*dudy*dvdx + dvdy*dvdy);
rhs = b*rho*deltaX*deltaX*deltaY*deltaY;
p_new(j,i) = ( (deltaY*deltaY)*(p(j,1)+p(j,i-1)) + (deltaX*deltaX)*(p(j+1,i)+p(j-1,i)) - rhs) / (2.0*(deltaX*deltaX + deltaY*deltaY));
i = 1;
dudx = (u(j,i+1)-u(j,nx))/(2*deltaX);
dvdx = (v(j,i+1)-v(j,nx))/(2*deltaX);
dudy = (u(j+1,i)-u(j-1,i))/(2*deltaY);
dvdy = (v(j+1,i)-v(j-1,i))/(2*deltaY);
b = ((1/deltaT)*(dudx + dvdy))-(dudx*dudx + 2*dudy*dvdx + dvdy*dvdy);
rhs = b*rho*deltaX*deltaX*deltaY*deltaY;
p_new(j,i) = ( (deltaY*deltaY)*(p(j,i+1)+p(j,nx)) + (deltaX*deltaX)*(p(j+1,i)+p(j-1,i)) - rhs) / (2.0*(deltaX*deltaX + deltaY*deltaY));
end
p_norm = abs(norm(p_new)-norm(p_old))/norm(abs(p_old));
count = count + 1;
end
end
Page updated
Google Sites
Report abuse