Incompressible Steady Flow(Couette Flow & Backward Facing Step) using MATLAB
In this project, Couette flow ( flow of a viscous fluid in the space between two surfaces, one of which is moving tangentially relative to the other) is being studied.The main objective of the project is to discretize the Navier-Stokes equation using Finite Difference method, and solve them numerically using SIMPLE (Semi-Implicit Method for Pressure Linked Equations), to analyse the flow of viscous fluid between two surfaces, in which the top surface is moving relative to the bottom plate.
The second part of the project is to analyse the flow over 2D Planar Backward Facing step using same Matlab code, only just changing the boundary conditions and initial conditions in the code.
Project Part1: Couette FlowIn fluid dynamics, Couette flow is defined as the flow of viscous fluid in the space between to surfaces, one of which is moving relative to the other, tangentially. The flow is driven by the viscous drag force acting on the fluid. This configuration can model certain practical problems, like flow in lightly loaded journal bearings, in viscometry and to demonstrate approximation of reversibility.
This type of flow is named in the honour of Maurice Couette, a Professor of Physics at French University of Angers in late 19th century.
Project Part 2: Backward Facing StepFor the second part of project, we will be solving fully developed laminar channel flow on a backward facing step as given in the below given figure:
Project Approach: SIMPLE method is employed for solving Navier-Stokes equation:
Staggered Grid setting:
Results and Discussion:
Project: Part1Starting from the assumed conditions for pressure and velocities at all points, variables are iterated to the correct values which satisfy the mass imbalance equation i.e. are actually physical.
In the below given curve, velocity component in x direction is plotted as function of vertical distance across the duct. Profiles are shown for various iteration numbers from 4 to 300 iteration steps.
As we can see in the above given graph, with increase in iteration step, values for x component of velocity-U reaches convergence to the actual physical value. In the Matlab code for the above given case, actual convergence occurs at iteration number-617 for which variation in maximum value of mass imbalance equation-‘d’<1e-5. Of course convergence reaches faster as the required maximum value for ‘ d’ is increased. For convergence around 300th iteration, value for ‘d’ comes out to be 4e-4.
Evolution of U with iterations:
Pressure contours for different iteration number:
Project Part 2: Backward Facing Step:
Results Discussion and Conclusion:Project Part1: Couette Flow
• Variation of velocity (x component of velocity) contours clearly illustrate the spread of velocity from top plate to bottom as the iterations are carried out till convergence(mass imbalance ‘d’-->0). With increase in iteration, velocity magnitude increases with each iteration near the middle region, as the effect of upper plate moving reaches towards the bottom plate.
• The velocity spike carried out at grid point (15,5) in the initial conditions, is visible at initial iterations in pressure plots, but as the iteration reaches convergence, because of the pressure correction method, that initial spike automatically reduces and spreads throughout the whole field and somewhere between 50 and 150 iterations, that spike completely disappears.
• Convergence occurs at iteration number-617 for which variation in maximum value of mass imbalance equation-‘d’<1e-5. For convergence around 300th iteration, value for ‘d’ comes out to be 4e-4. Number of iterations also depends upon the convergence criteria set for pressure correction i.e. solving poisson equation for pressure.
• At around 300th iteration, the maximum value for d=4.2e-4, but for convergence d<1e-5, occurs after 617 iterations, with max d=9e-6. Similarly for 150 iterations, max d= 1e-3. Hence it depends upon the application, the kind of accuracy is required. Lower the value of d, higher are the number of iterations required.
Project Part2: Backward Facing step
• Velocity contours(x component of velocity-U) clearly illustrates the variation of velocity at each point in domain as iteration reaches convergence. The hydraulic domain at the lower part of inflow shows negative magnitude , that means the direction of flow reverses at the lower part, below inflow which illustrates the effect of backward facing step. This can also be verified in the velocity vector diagram.
• With increase in the number of grid points , the value of dt had to be decreased for stable convergence, because the grid spacing decreases as the number of points are increased.
• Number of iterations increased as the Reynold number is decreased for the same mass imbalance(dmax=1e-5). For instance in case of grid(51X21), for Reynold number-1000, number of iterations for main loop is 89, but for Reynold number-10, number of iterations shoots to about 675 iterations.
• With increase in the number of grid points(101X41), value of dmax had to be increased, as the it takes much more time to convergence for the same dmax value as that for smaller grid(51X21).
clear all
close all
clc
%%%%%given%%%%%%%%%
dt=1e-3;
dens=0.002377; %slug/ft^2 ---incompressible flow so always constant
mu=(dens*0.01*1)/63.6;
%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%Grid%%%%%%%
D=[0.5 0.01]; %Domain size in x and y direction
Xp=21; % number of Pressure points in x direction
Yp=11; % number of Pressure points in y direction
Xu=Xp+1; % number of points for U-velocity in x direction
Yu=Yp; % number of points for U-velocity in y direction
Xv=Xp+2; %number of points for V-velocity in x direction
Yv=Yp+1; %number of points for V-velocity in y direction
dx=D(1)/(Xp-1); %grid spacing in x direction
dy=D(2)/(Yp-1);%grid spacing in y direction
y_len=linspace(0,0.01,11);
P=zeros(Yp,Xp);%array formation
U=zeros(Yu,Xu);%array formation
V=zeros(Yv,Xv);%array formation
d=zeros(Yp,Xp);%array formation
P_star=zeros(Yp,Xp); %initial pressure guess for iteration
P_dash=zeros(Yp,Xp); % P' for pressure poisson equation
%Boundary conditions
%%%%%% wall %%%%%%%
for i=1:Xu
U(1,i)=1; % U-velocity at upper wall
V(1,i)=0; % V-velocity at upper wall
U(Yu,i)=0; % U-velocity at bottom wall
V(1,i)=0; % V-velocity at bottom wall
end
%%%%%%%%%%%%%%%%%%%
%%%%%%% inflow %%%%%%%%%
for i=1:Yp
P_dash(i,1)=0; %P' at inflow boundary
end
for i=1:Yv
V(i,2)=0; % V at inflow boundary V(i,1)---ghost point
end
%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% outflow %%%%%%%%
for i=1:Yp
P_dash(i,Xp)=0; %P' at inflow boundary
end
%%%%%%%%%%%%%%%%%%%%%%%%
%%% initial velocity spike for V %%%
V(5,15)=0.5; % initial velocity spike
%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Density-velocity matrix %%%%%%%%%
rhou=dens*U;
rhov=dens*V;
rho_u2=dens*(U.*U);
rho_v2=dens*(V.*V);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
P_dash_new=P_dash;
%%%%%%%%%%%%%for checking-code%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% Momentum Loop for u* v* %%%%%%%%%%
% A_s=-(((rho_u2(2,3)-rho_u2(2,1))/(2*dx))+(((rhou(1,2)*(V(2,2)+V(2,3))*0.5)-(rhou(3,2)*(V(3,2)+V(3,3))*0.5))/(2*dy)))+mu*(((U(2,3)-2*U(2,2)+U(2,1))/(dx)^2)+((U(1,2)-2*U(2,2)+U(3,2))/(dy)^2));
%
% rhou(2,2)=rhou(2,2)+A_s*dt-(dt/dx)*(P(2,2)-P(2,1));
% U(2,2)=rhou(2,2)/dens;
% rho_u2(2,2)=rhou(2,2)*U(2,2);
%
% A_s1=-(((rho_u2(2,4)-rho_u2(2,2))/(2*dx))+(((rhou(1,3)*(V(2,3)+V(2,4))*0.5)-(rhou(3,3)*(V(3,3)+V(3,4))*0.5))/(2*dy)))+mu*(((U(2,4)-2*U(2,3)+U(2,2))/(dx)^2)+((U(1,3)-2*U(2,3)+U(3,3))/(dy)^2));
%
% rhou(2,3)=rhou(2,3)+A_s1*dt-(dt/dx)*(P(2,3)-P(2,2));
%%%-- x---i+1=2; y----j+1=2;,,in matlab its inverse of cartesian coordinates
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%5for plotting u vs iteration no.%%%%%%%%%%%%%%%
% count=0;
% for t = 1:5
%
% count = count+1;
%
% l = [4 20 50 150 300];
%
% for iter = 1 : l(t)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dmax=1;itrd=0; %for iterations
%%%%main loop%%%%%
while(dmax>1e-5)
%%%%velocity loop%%%
for j=1:Yp-2
for i=1:Xp-1
A_s(j+1,i+1)=-(((rho_u2(j+1,i+2)-rho_u2(j+1,i))/(2*dx))+(((rhou(j,i+1)*(V(j+1,i+1)+V(j+1,i+2))*0.5)-(rhou(j+2,i+1)*(V(j+2,i+1)+V(j+2,i+2))*0.5))/(2*dy)))+mu*(((U(j+1,i+2)-2*U(j+1,i+1)+U(j+1,i))/(dx)^2)+((U(j,i+1)-2*U(j+1,i+1)+U(j+2,i+1))/(dy)^2));
rhou(j+1,i+1)=rhou(j+1,i+1)+A_s(j+1,i+1)*dt-(dt/dx)*(P_star(j+1,i+1)-P_star(j+1,i));
U(j+1,i+1)=rhou(j+1,i+1)/dens;
rho_u2(j+1,i+1)=rhou(j+1,i+1).*U(j+1,i+1);
end
end
U(:,1)=U(:,2);
U(:,Xu)=U(:,Xu-1);
for j=1:Yp-1
for i=1:Xp
B_s(j+1,i+1)=-(((rho_v2(j,i+1)-rho_v2(j+2,i+1))/(2*dy))+(((rhov(j+1,i+2)*(U(j+1,i+1)+U(j,i+1))*0.5)-(rhov(j+1,i)*(U(j,i)+U(j+1,i))*0.5))/(2*dx)))+mu*(((V(j,i+1)-2*V(j+1,i+1)+V(j+2,i+1))/(dy)^2)+((V(j+1,i+2)-2*V(j+1,i+1)+V(j+1,i))/(dx)^2));
rhov(j+1,i+1)=rhov(j+1,i+1)+(B_s(j+1,i+1)*dt)-(dt/dy)*(P_star(j,i)-P_star(j+1,i));
V(j+1,i+1)=rhov(j+1,i+1)/dens;
rho_v2(j+1,i+1)=rhov(j+1,i+1).*V(j+1,i+1);
end
end
%%%%%velocity loop%%%%%%%
%%%%for checking%%%%%%%
% %for last line of points for V---in x direction
% for i=1:Xp
% B_s=-(((rho_v2(Yp-1,i+1)-rho_v2(Yp+1,i+1))/(2*dy))+(((rhov(Yp,i+2)*(U(Yp,i+1)+U(Yp-1,i+1))*0.5)-(rhov(Yp,i)*(U(Yp-1,i)+U(Yp,i))*0.5))/(2*dx)))+mu*(((V(Yp-1,i+1)-2*V(Yp,i+1)+V(Yp+1,i+1))/(dy)^2)+((V(Yp,i+2)-2*V(Yp,i+1)+V(Yp,i))/(dx)^2));
% rhov(Yp,i+1)=rhov(Yp,i+1)+B_s*dt-(dt/dy)*(P_star(Yp-1,i)-P_star(Yp,i));
% V(Yp,i+1)=rhov(Yp,i+1)/dens;
% rho_v2(Yp,i+1)=rhov(Yp,i+1).*V(Yp,i+1);
% end
%
% % for last line of points for V---in y direction
% for j=1:Yp-1
% B_s=-(((rho_v2(j,Xp+1)-rho_v2(j+2,Xp+1))/(2*dy))+(((rhov(j+1,Xp+2)*(U(j+1,Xp+1)+U(j,Xp+1))*0.5)-(rhov(j+1,Xp)*(U(j,Xp)+U(j+1,Xp))*0.5))/(2*dx)))+mu*(((V(j,Xp+1)-2*V(j+1,Xp+1)+V(j+2,Xp+1))/(dy)^2)+((V(j+1,Xp+2)-2*V(j+1,Xp+1)+V(j+1,Xp))/(dx)^2));
% rhov(j+1,Xp+1)=rhov(j+1,Xp+1)+B_s*dt-(dt/dy)*(P_star(j,Xp)-P_star(j+1,Xp));
% V(j+1,Xp+1)=rhov(j+1,Xp+1)/dens;
% rho_v2(j+1,Xp+1)=rhov(j+1,Xp+1).*V(j+1,Xp+1);
% end
%%%%%%%%%%%%%%%%%%%%
% modification values for velocities at boundaries
% U(:,1)=U(:,2);
% U(:,Xu)=U(:,Xu-1);
%V(:,1)=V(:,2);
V(:,Xv)=V(:,Xv-1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i= 2:Xp-1
for j=2:Yp-1
d(j,i)=(1/dx)*(rhou(j,i+1)-rhou(j,i))+(1/dy)*(rhov(j,i+1)-rhov(j+1,i+1));
end
end
a=2*dt*(((1/dx)^2)+((1/dy)^2));
b=-dt/(dx)^2;
c=-dt/(dy)^2;
%%%%%%%% Poisson Equation-Loop for Pressure %%%%%%%%%%%%%
error=1;
itr=0;
while(error>1e-10)
for i= 2:Xp-1
for j=2:Yp-1
P_dash_new(j,i)=(-1/a)*(b*P_dash(j,i+1)+b*P_dash_new(j,i-1)+c*P_dash(j+1,i)+c*P_dash_new(j-1,i)+d(j,i));
E(j,i)=abs(P_dash_new(j,i)-P_dash(j,i));
end
end
error=max(E(:));
%P_dash=P_dash+1.1*(P_dash_new-P_dash);
P_dash=P_dash_new;
itr=itr+1;
end
P_star=P_star+0.1*P_dash;
d1=abs(d);
dmax=max(d1(:));
%%%%%plotting -Pressure and velocity at different iteration points
itrd=itrd+1;
if(itrd>500)
break;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% pcolor(flipud(P_star));shading interp ;
plot(flipud(P_star(:,19)),y_len)
pause(0.1);
axis([0 1e-5 0 0.01])
%colorbar
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% u_plot(:,count) = flipud(U(:,15));
%
% end
% figure();
% %
% plot((u_plot(:,1)), y_len','--*r','Markerface','r');
%
% hold on;
%
% plot((u_plot(:,2)), y_len','--*r');
%
% plot((u_plot(:,3)), y_len','--*b');
%
% plot((u_plot(:,4)), y_len','--*g');
%
% plot((u_plot(:,5)), y_len','--ok');
%
% %set(gca,'fontsize',12,'fontweight','bold')
%
% legend('k = 4','k = 20','k = 50','k = 150','k = 300')
%
% xlabel('u, ft/s');
%
% ylabel('y, ft');
%
% %axis([0 1 0 H]);
% grid on;
% % % %Grid plotting
% %
% % [gxp gyp] = meshgrid(0:dx:0.5,0:dy:0.01);
% % [gxu gyu] = meshgrid(-0.5*dx:dx:0.5+0.5*dx,0:dy:0.01);
% % [gxv gyv] = meshgrid(-1*dx:dx:0.5+1*dx,-0.5*dy:dy:0.01+0.5*dy);
% %
% % plot(gxp,gyp,'sk',gxu,gyu,'or',gxv,gyv,'xg')
% %
%Project case 2:
clear all
close all
clc
%given
dt=1e-3; %e-3 %with more points dt value decreases for convergence
dens=0.002377; %slug/ft^2 ---incompressible flow so always constant
mu=(dens*2*1)/1000; %%%----check %%63.6--10;1000; 10--675iteration; 1000-89
%%%
D=[4.5 2]; %Domain size in x and y direction
Xp=51; % number of Pressure points in x direction %51
Yp=21; % number of Pressure points in y direction %21
Xu=Xp+1; % number of points for U-velocity in x direction
Yu=Yp; % number of points for U-velocity in y direction
Xv=Xp+2; %number of points for V-velocity in x direction
Yv=Yp+1; %number of points for V-velocity in y direction
dx=D(1)/(Xp-1);
dy=D(2)/(Yp-1);
y_len=linspace(0,2,Yu);
x_len=linspace(0,4.5,Xu);
P=zeros(Yp,Xp);
U=zeros(Yu,Xu);
V=zeros(Yv,Xv);
d=zeros(Yp,Xp);
P_star=zeros(Yp,Xp); %initial pressure guess for iteration
P_dash=zeros(Yp,Xp); % P' for pressure poisson equation
%Boundary conditions
%%%%%% wall %%%%%%%
for i=1:Xu
U(1,i)=0; % U-velocity at upper wall---make this 0 for optional case
V(1,i)=0; % V-velocity at upper wall-----check
U(Yu,i)=0; % U-velocity at bottom wall
V(1,i)=0; % V-velocity at bottom wall-----check
end
%%%%%%%%%%%%%%%%%%%
itr=1;
%%%%%%% inflow %%%%%%%%%--- this we have to change for U
for i=1:Yp
P_dash(i,1)=0; %P' at inflow boundary
end
r=linspace(1,0,((Yu-1)/4)+1);
r1=linspace(0,1,((Yu-1)/4)+1);
for i= 1:((Yu-1)/4)+1%for inflow -middle of inlet
U(i,1)=1*(1-(r(i))^2);
end
for i= ((Yu-1)/4)+2:((Yu-1)/2)+1%for inflow -middle of inlet to end middle
U(i,1)=U(i-2*itr,1);
itr=itr+1;
end
for i=1:Yv
V(i,2)=0; % V at inflow boundary V(i,1)---ghost point
end
%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% outflow %%%%%%%%
for i=1:Yp
P_dash(i,Xp)=0; %P' at inflow boundary
end
%%%%%%%%%%%%%%%%%%%%%%%%
%%% initial velocity spike for V %%%
% V(5,15)=0.5; % initial velocity spike
%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Density-velocity matrix %%%%%%%%%
rhou=dens*U;
rhov=dens*V;
rho_u2=dens*(U.*U);
rho_v2=dens*(V.*V);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
P_dash_new=P_dash;
dmax=1;itrd=0;
%for p=1:300;
while(dmax>1e-5)
for j=1:Yp-2
for i=1:Xp-1
A_s(j+1,i+1)=-(((rho_u2(j+1,i+2)-rho_u2(j+1,i))/(2*dx))+(((rhou(j,i+1)*(V(j+1,i+1)+V(j+1,i+2))*0.5)-(rhou(j+2,i+1)*(V(j+2,i+1)+V(j+2,i+2))*0.5))/(2*dy)))+mu*(((U(j+1,i+2)-2*U(j+1,i+1)+U(j+1,i))/(dx)^2)+((U(j,i+1)-2*U(j+1,i+1)+U(j+2,i+1))/(dy)^2));
rhou(j+1,i+1)=rhou(j+1,i+1)+A_s(j+1,i+1)*dt-(dt/dx)*(P_star(j+1,i+1)-P_star(j+1,i));
U(j+1,i+1)=rhou(j+1,i+1)/dens;
rho_u2(j+1,i+1)=rhou(j+1,i+1).*U(j+1,i+1);
end
end
U(:,Xu)=U(:,Xu-1);
for j=1:Yp-1
for i=1:Xp
B_s(j+1,i+1)=-(((rho_v2(j,i+1)-rho_v2(j+2,i+1))/(2*dy))+(((rhov(j+1,i+2)*(U(j+1,i+1)+U(j,i+1))*0.5)-(rhov(j+1,i)*(U(j,i)+U(j+1,i))*0.5))/(2*dx)))+mu*(((V(j,i+1)-2*V(j+1,i+1)+V(j+2,i+1))/(dy)^2)+((V(j+1,i+2)-2*V(j+1,i+1)+V(j+1,i))/(dx)^2));
rhov(j+1,i+1)=rhov(j+1,i+1)+(B_s(j+1,i+1)*dt)-(dt/dy)*(P_star(j,i)-P_star(j+1,i));
V(j+1,i+1)=rhov(j+1,i+1)/dens;
rho_v2(j+1,i+1)=rhov(j+1,i+1).*V(j+1,i+1);
end
end
V(:,Xv)=V(:,Xv-1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i= 2:Xp-1
for j=2:Yp-1
d(j,i)=(1/dx)*(rhou(j,i+1)-rhou(j,i))+(1/dy)*(rhov(j,i+1)-rhov(j+1,i+1));
end
end
a=2*dt*(((1/dx)^2)+((1/dy)^2));
b=-dt/(dx)^2;
c=-dt/(dy)^2;
%%%%%%%% Poisson Equation for Pressure %%%%%%%%%%%%%
error=1;
itr=0;
while(error>1e-10)
for i= 2:Xp-1
for j=2:Yp-1
P_dash_new(j,i)=(-1/a)*(b*P_dash(j,i+1)+b*P_dash_new(j,i-1)+c*P_dash(j+1,i)+c*P_dash_new(j-1,i)+d(j,i));
E(j,i)=abs(P_dash_new(j,i)-P_dash(j,i));
end
end
error=max(E(:));
%P_dash=P_dash+1.1*(P_dash_new-P_dash);
P_dash=P_dash_new;
itr=itr+1;
end
P_star=P_star+0.1*P_dash;
d1=abs(d);
dmax=max(d1(:));
itrd=itrd+1;
% if(itrd>300)
% break;
% end
%plot(flip(U(:,3)),y_len')
%pcolor((flipud(U(:,2:Xu))));shading interp ;
%pcolor((flipud(V(2:21,2:51))));shading interp ;
%pcolor((flipud(P_star)));shading interp ;
%pause(0.001)
%colorbar
%axis([-0.5 1 0 2])
% hold on
end
%plot(flip(P_star(:,14)),y_len')
%
%%%%%%%%%%%Vector plots%%%%%%%%%%%%
y_len=linspace(0,2,Yp);
x_len=linspace(0,4.5,Xp);
[Xg,Yg]=meshgrid(x_len,y_len);
quiver(Xg,flip(Yg),U(:,2:Xu),V(2:Yv,2:Xv-1),10)
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% axis([-1 4.5 0 2.5])