Backward Facing Step using Simple Algorithm
 

Hi all,

I've solved the lid-driven cavity flow problem in MATLAB using the Simple algorithm. I am now working to solve the backward-facing step flow problem, also using the simple algorithm. Based on my literature review, I should be able to modify the lid-driven cavity flow code for this purpose by simply changing the boundary conditions. I've tried to do this several different ways, but can't seem to get it right.

Here are the boundary conditions:

inlet:
u = 4*Umax*(y/Hs)*(1-y/Hs); Hs = height of parabolic flow inlet
v = 0

outlet:
du/dx = dv/dx = 0

top and bottom walls:
u = v = 0

Reynolds number is defined as Re = Umax(H-Hs)/v; H = total height

I've attached the MATLAB code that I've been working on. Any help would be greatly appreciated!

It is much better you explain the problems and post some plot of your solution

Okay, here are some of the simulation details:

H = 1.5
Hs = 1.0
h = H-Hs
L = 30*h
Re = 10 (uMax = 20)
Iterations = 1000

Plot is attached. Notice that all of the velocities are negative, which is counter-intuitive with a positive x parabolic inflow. As I increase the Re to 100, the u-velocity becomes more negative.

since H=1.5 and Hs=1, using you function the flow inlet u(y) is described for y=1 to 1.5 and u is negative at the inflow

Thanks FMDenaro.

I was given that specific parabolic function as part of the problem description, and had noticed that the velocity profile was negative at the left boundary. I think I will either have to change the origin point from the bottom left to top left, or derive a new inlet parabolic function.

Okay, next question:

To solve for x-momentum as uStar, I use the following code:

% setup coefficients
for i = 2:Nx
for j = 2:Ny+1
% convective flux using central differencing
uCe = dx*0.5*(uOld(i,j)+uOld(i+1,j));
uCw = dx*0.5*(uOld(i,j)+uOld(i-1,j));
uCn = dy*0.5*(vOld(i,j)+vOld(i+1,j));
uCs = dy*0.5*(vOld(i,j-1)+vOld(i+1,j-1));
% central coefficient using applied under-relaxation
uDx = dy/dx*(1/Re);
uDy = dx/dy*(1/Re);
auP(i,j) = (0.5*(uCe+uCw+uCn+uCs)+2*uDx+2*uDy)/alphaU;
% boundary coefficients using hybrid schemes
% change numerical method as required
auE_h = [-uCe,-0.5*uCe+uDx,0];
auW_h = [uCw,-0.5*uCw+uDx,0];
auN_h = [-uCn,-0.5*uCn+uDy,0];
auS_h = [uCs,-0.5*uCs+uDy,0];
auE(i,j) = max(auE_h);
auW(i,j) = max(auW_h);
auN(i,j) = max(auN_h);
auS(i,j) = max(auS_h);
% set u velocity component of PCE
dU(i,j) = dx/auP(i,j);
end
end
% use previous calculation as initial guess in numerical method
uStar = uOld;

(similar code is used to solve for y-momentum).

Notice the uDx and uDy values contain the Re value. When I have it coded as is, the solution blows up, as in the first image attached. However, if I do not divide by Re (i.e. for uDx = dy/dx*Re), then my solution trends in the right way, but with the wrong u-momentum values (as in the second image attached).

Note that for the solution that 'blows up' I only ran 20 iterations because it fails with more, while for the solution that trends properly with wrong values, I ran at 200 and 1000 iterations (i.e. it doesn't 'blow up' with more iterations). I've also re-attached my code, as it has been changed since the first post.

Note also that when coded as uDx = dy/dx*(1/Re), the lid-driven cavity problem solves just fine.

Again, any advice would be appreciated.

I've continued to tinker with the backward step flow MATLAB code, and the simulation is trending in the right way.

My question now is about boundary conditions. I have the following:

Inlet:
u = parabolic
v = 0

Outlet:
du/dx = dv/dx = 0

Top and Bottom:
u = v = 0

I've coded this as follows:

%Left and Right
for j = 1:Ny
if y(1,j)>Hs&&y(1,j)<H
Left(1,j) = 4*uMax*(y(1,j)/Hs)*(1-(y(1,j)/Hs));
%Left(1,j) = 100;
else
Left(1,j) = 0.0;
end
uOld(1,j) = uOld(2,j) - Left(1,j);
uOld(Nx+2,j) = uOld(Nx+1,j);
%vOld(1,j) = vOld(2,j);
%vOld(Nx+2,j) = vOld(Nx+1,j);
end

%Top and Bottom
for i = 1:Nx+1
%uOld(i,1)= uOld(i,2);
%uOld(i,Ny+1)= uOld(i,Ny);
vOld(i,1)= vOld(i,2);
vOld(i,Ny+1)= vOld(i,Ny);
end

Is this a reasonable approach?

I don't think that posting your code is useful... if you have a problem in the solution describe it, give details, figures and a list of the main setting.

Thanks FMDenaro. I don't intend to use CFD Online as a 'code editor' - I just assumed that someone may have done something similar to me and could offer quick advice.

My assumption is that I do not need to change the SIMPLE algorithm code from that used in the lid-driven problem, and instead need to change my left and right (inlet and outlet) boundary conditions for the backward facing step problem.

My approach (given the boundary conditions presented in my previous post) has been to define the parabolic function at x = 1 for all y using a 'for' statement. I've then attempted to code the boundary conditions as previously described.

My results are okay, and I can 'tune' them using the the kinematic viscosity. However, my concern is that my original assumption (i.e. that the SIMPLE algorithm doesn't change) is incorrect.

If you wrote in a general way, the algorithm incorporates general BC.s so that you have only to set BC.s for inflow+wall on the left side and an outflow on the right side.
What I suggest is to debug the code, sometimes a code can work fine for totally confined flows but can show problems for inflow/outflow conditions. And, often, a bug is hidden...

Very helpful FMDenaro - thanks! I've actually noticed several bugs in the confined flow scenario that needed fixing for the inflow/outflow scenario.

I'll keep working on the problem and post when I have more questions or a consistent solution.

I've finished working on this bit of research, and I've posted the resulting paper on Scribd:

https://www.scribd.com/document/3563...MPLE-Algorithm

I'll be converting the MATLAB code here developed to Python code in the upcoming months, so I'll most likely post additional updates. Stay tuned!

Great
 

Can you send me
How do you write this code
I want your file
Mohammadaminjafari1372@yahoo.com

Thank you for sharing the paper.

Dear Aaron: thanks for posting the code..

minor comment on your code:



should your (yy) in the code be :
yy= H-dy/2:-dy:dy/2
instead of

yy= H:-dy:0

also the same change for y ?