
[Sponsors] 
March 6, 2017, 12:12 
Backward Facing Step using Simple Algorithm

#1 
New Member
aaron outhwaite
Join Date: Mar 2016
Posts: 10
Rep Power: 6 
Hi all,
I've solved the liddriven cavity flow problem in MATLAB using the Simple algorithm. I am now working to solve the backwardfacing step flow problem, also using the simple algorithm. Based on my literature review, I should be able to modify the liddriven 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)*(1y/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(HHs)/v; H = total height I've attached the MATLAB code that I've been working on. Any help would be greatly appreciated! 

March 6, 2017, 12:17 

#2 
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 5,342
Rep Power: 57 
It is much better you explain the problems and post some plot of your solution


March 6, 2017, 12:31 

#3 
New Member
aaron outhwaite
Join Date: Mar 2016
Posts: 10
Rep Power: 6 
Okay, here are some of the simulation details:
H = 1.5 Hs = 1.0 h = HHs L = 30*h Re = 10 (uMax = 20) Iterations = 1000 Plot is attached. Notice that all of the velocities are negative, which is counterintuitive with a positive x parabolic inflow. As I increase the Re to 100, the uvelocity becomes more negative. 

March 6, 2017, 12:38 

#4 
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 5,342
Rep Power: 57 
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: 6 
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: 6 
Okay, next question:
To solve for xmomentum 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(i1,j)); uCn = dy*0.5*(vOld(i,j)+vOld(i+1,j)); uCs = dy*0.5*(vOld(i,j1)+vOld(i+1,j1)); % central coefficient using applied underrelaxation 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 ymomentum). 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 umomentum 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 reattached my code, as it has been changed since the first post. Note also that when coded as uDx = dy/dx*(1/Re), the liddriven cavity problem solves just fine. Again, any advice would be appreciated. 

March 9, 2017, 10:04 

#7 
New Member
aaron outhwaite
Join Date: Mar 2016
Posts: 10
Rep Power: 6 
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? 

March 9, 2017, 11:40 

#8 
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 5,342
Rep Power: 57 
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: 6 
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 liddriven 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: 5,342
Rep Power: 57 
Quote:
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: 6 
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: 6 
I've finished working on this bit of research, and I've posted the resulting paper on Scribd:
https://www.scribd.com/document/3563...MPLEAlgorithm 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! 

November 28, 2017, 15:53 
Great

#13 
New Member
Ebara
Join Date: Nov 2017
Posts: 1
Rep Power: 0 

March 10, 2018, 12:43 

#14  
New Member
Darsh Nathawani
Join Date: Feb 2018
Posts: 4
Rep Power: 4 
Quote:


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= Hdy/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 
Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Backward facing step; x/H; U/U0; diemsionless analysis; paraview  idrama  OpenFOAM  1  December 3, 2013 17:10 
Micro Scale Pore, icoFoam  gooya_kabir  OpenFOAM Running, Solving & CFD  2  November 2, 2013 13:58 
Velocity drop over backward facing step  warlocklw  CFX  2  February 9, 2012 21:24 
airfuel mixing in backward facing step  abdullahkarimi  Main CFD Forum  0  March 19, 2011 18:39 
Corner Vortex for Backward Facing Step  Patrick G. Hu  Main CFD Forum  0  August 6, 2002 09:34 