# Backward Facing Step using Simple Algorithm

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

March 6, 2017, 12:12
Backward Facing Step using Simple Algorithm
#1
New Member

aaron outhwaite
Join Date: Mar 2016
Posts: 10
Rep Power: 10
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!
Attached Files
 Backward_AO_CFD_Online.zip (2.2 KB, 139 views)

 March 6, 2017, 12:17 #2 Senior Member   Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,777 Rep Power: 71 It is much better you explain the problems and post some plot of your solution lcarasik and TheVictorVS like this.

March 6, 2017, 12:31
#3
New Member

aaron outhwaite
Join Date: Mar 2016
Posts: 10
Rep Power: 10
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.
Attached Images
 u_centerline_Re_10.jpg (15.1 KB, 69 views)

 March 6, 2017, 12:38 #4 Senior Member   Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,777 Rep Power: 71 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

 March 6, 2017, 12:47 #5 New Member   aaron outhwaite Join Date: Mar 2016 Posts: 10 Rep Power: 10 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.

March 6, 2017, 14:54
#6
New Member

aaron outhwaite
Join Date: Mar 2016
Posts: 10
Rep Power: 10
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.
Attached Images
 u_centerline_Re_100_fail.jpg (15.9 KB, 34 views) u_centerline_Re_100_fail_trend.jpg (14.8 KB, 35 views)
Attached Files
 Backward_Step_Matlab.zip (2.2 KB, 65 views)

 March 9, 2017, 10:04 #7 New Member   aaron outhwaite Join Date: Mar 2016 Posts: 10 Rep Power: 10 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)

 March 9, 2017, 11:40 #8 Senior Member   Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,777 Rep Power: 71 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.

 March 9, 2017, 12:49 #9 New Member   aaron outhwaite Join Date: Mar 2016 Posts: 10 Rep Power: 10 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.

March 9, 2017, 12:55
#10
Senior Member

Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,777
Rep Power: 71
Quote:
 Originally Posted by aar_out 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...

 March 9, 2017, 13:07 #11 New Member   aaron outhwaite Join Date: Mar 2016 Posts: 10 Rep Power: 10 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.

 August 15, 2017, 19:41 #12 New Member   aaron outhwaite Join Date: Mar 2016 Posts: 10 Rep Power: 10 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! lcarasik, jhawisa, MMNCH and 1 others like this.

 November 28, 2017, 15:53 Great #13 New Member   Ebara Join Date: Nov 2017 Posts: 1 Rep Power: 0 Can you send me How do you write this code I want your file Mohammadaminjafari1372@yahoo.com

March 10, 2018, 12:43
#14
New Member

Darsh Nathawani
Join Date: Feb 2018
Posts: 4
Rep Power: 8
Quote:
 Originally Posted by aar_out 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!
Thank you for sharing the paper.

 March 25, 2019, 16:15 #15 New Member   Jamal Join Date: Jul 2017 Posts: 2 Rep Power: 0 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 ?

 Tags backward facing step, boundary condition, matlab, simple algorithm