# Matlab code for pipe flow

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

 August 31, 2011, 13:24 #21 New Member   Join Date: Jul 2011 Posts: 25 Rep Power: 8 Sponsored Links Hi there, I defined the boundary conditions for the outlet like this: P = 0, du/dn=0, dv/dn=0 As to the error I mentioned, I still think I am right. Regards

 September 1, 2011, 00:18 #22 Member   Prasanth P Join Date: May 2009 Posts: 30 Rep Power: 10 Hi, I went through the code. I think you are right regarding the error. I guess the issue is with the Poisson solver that we are using. I am using Gauss JAcobi method. What method are you using for solving the Poisson equation?

 September 5, 2011, 13:12 #24 New Member   Join Date: Jul 2011 Posts: 25 Rep Power: 8 Hi Vincent, I think I can answer your questions about U+Ubc. If you work through the code, you will find that Ubc is not actually for applying the boundary conditions in the normal sense, it is really for computing the RHS, or say the source term of the system of equations, which is U+Ubc. Notice that in the notes about this code, the second step is for implicit viscosity, whose formulation is (U** - U*)/dt = (Uxx** + Uyy**)/Re, simple rearrangement of terms yields U** - dt/Re*(Uxx** + Uyy**) = U*. The U* is the U in U+Ubc. How about Ubc? If we implement the above formulation to the boundary nodes, some terms need to be moved to the RHS too since some unknowns for interior nodes become known for boundary nodes. Therefore the final complete RHS for solving this ystem is U+Ubc. However, the error I described in previous thread does exist. Hope this helps. If not clear enough, I can type in the complete formulation if I've got time. Cici

 September 6, 2011, 00:10 Regarding the code #25 Member   Prasanth P Join Date: May 2009 Posts: 30 Rep Power: 10 Hi Vincent, "Secondly I had already noticed that the code only works correctly for hx=hy. Though this might not be been intended, it is faster to calculate on a completely regular grid. So this potential bug does not have to cause any problems while using the code." This bug is because of the Ubc definition. Try writing the Ubc as posted by Cici you would find that the code would start working for hx/=hy. This bug arises due to the mixup of definitions of Ubc and Vbc. Cheers PP

September 6, 2011, 03:19
#26
New Member

Vincent
Join Date: Jul 2011
Posts: 29
Rep Power: 8
Quote:
 Originally Posted by prasanthnitt Hi Vincent, "Secondly I had already noticed that the code only works correctly for hx=hy. Though this might not be been intended, it is faster to calculate on a completely regular grid. So this potential bug does not have to cause any problems while using the code." This bug is because of the Ubc definition. Try writing the Ubc as posted by Cici you would find that the code would start working for hx/=hy. This bug arises due to the mixup of definitions of Ubc and Vbc. Cheers PP
Well, that is a really clear confirmation of the error made defining Ubc and Vbc. Thanks for figuring that one out.

Quote:
 Hi Vincent, I think I can answer your questions about U+Ubc. If you work through the code, you will find that Ubc is not actually for applying the boundary conditions in the normal sense, it is really for computing the RHS, or say the source term of the system of equations, which is U+Ubc. Notice that in the notes about this code, the second step is for implicit viscosity, whose formulation is (U** - U*)/dt = (Uxx** + Uyy**)/Re, simple rearrangement of terms yields U** - dt/Re*(Uxx** + Uyy**) = U*. The U* is the U in U+Ubc.
This part was completely clear, I already formulated this in my own post ("the matrix multiplication is (I +R(K2D))U^{**}=U^{*}).

Quote:
 How about Ubc? If we implement the above formulation to the boundary nodes, some terms need to be moved to the RHS too since some unknowns for interior nodes become known for boundary nodes. Therefore the final complete RHS for solving this ystem is U+Ubc. However, the error I described in previous thread does exist. Hope this helps. If not clear enough, I can type in the complete formulation if I've got time. Cici
I can't say I understand it yet, if you could derive why Ubc has to have
that expression I would be very thankfull.

U** -dt/Re(Uxx**+Uyy**)=U*
U** -dt/Re((Ui-1,j**-2Ui,j** +Ui+1,j**)/hx^2+(Ui,j-1** -2Ui,j**+Ui,j+1**)/hy^2)=U*

Now I understand that Ui-1,j, Ui+1,j, Ui,j-1 and Ui,j+1 are not defined at every pont. That is why we move some points the the RHS of this equation.

So for instance we add dt/Re(Ui-1,j)/hx^2 to the left side of our matrix U(1,:). If we correct the mistake I understand that this is the same as uW, ditto for uE.

Why, however, it is in the case of hy^2 terms that a 2 appears?
Thanks for the help,

Regards,

Vincent

PS: You can improve the code by rewriting the function B:

function B = avg(A,k)
if nargin<2, k = 1; end
if
size(A,1)==1, A = A'; end
if
k<2, B = conv2(A,[1;1]/2,'valid'); else B = avg(A,k-1); end
if
size(A,2)==1, B = B'; end

You will see that this is a faster way of calculating the same thing (makes quite a difference for large matrices).

PS2: wouldn't it be nice if this forum had a latex compiler?

 September 6, 2011, 05:16 Code and Pressure Bcn #27 Member   Prasanth P Join Date: May 2009 Posts: 30 Rep Power: 10 Hello Vincent, The correct expression would be: Ubc = dt/Re*([2*uS(2:end-1)' zeros(nx-1,ny-2) 2*uN(2:end-1)']/hy^2+... [uW;zeros(nx-3,ny);uE]/hx^2); Vbc = dt/Re*([vS' zeros(nx,ny-3) vN']/hy^2+... [2*vW(2:end-1);zeros(nx-2,ny-1);2*vE(2:end-1)]/hx^2); A 2 appears not only in the hy^2 terms, it also appears in hx^2 terms. For Ubc its in hy^2 term but in Vbc it is in hx^2 terms(from the above expressions). It appears so because in x direction the x velocity is staggered from the pressure variable and in the y direction, the y velocity is staggered from pressure. I hope this clears your doubt. Now getting back to the issue with the boundary conditions for pressure. Have any of you written an iterative solver for the pressure poisson equation? Dr. Seibold's code uses a direct method to solve the the poisson equation for the pressure. I am facing issue with my iterative procedure for solving the equations. Plus, when the out flow boundary condition is partial(u)/partial(n) = 0, how can you say that P=0 or dp/dn=0 is a correct boundary condition? Cheers PP

September 6, 2011, 05:56
#28
New Member

Vincent
Join Date: Jul 2011
Posts: 29
Rep Power: 8
Quote:
 Hello Vincent, The correct expression would be: Ubc = dt/Re*([2*uS(2:end-1)' zeros(nx-1,ny-2) 2*uN(2:end-1)']/hy^2+... [uW;zeros(nx-3,ny);uE]/hx^2); Vbc = dt/Re*([vS' zeros(nx,ny-3) vN']/hy^2+... [2*vW(2:end-1);zeros(nx-2,ny-1);2*vE(2:end-1)]/hx^2); A 2 appears not only in the hy^2 terms, it also appears in hx^2 terms. For Ubc its in hy^2 term but in Vbc it is in hx^2 terms(from the above expressions). It appears so because in x direction the x velocity is staggered from the pressure variable and in the y direction, the y velocity is staggered from pressure. I hope this clears your doubt.
Alright I understand that we defined U and V velocity fields on a staggered grid. In the nonstaggered direction we do not have a factor, as we would expect.
In the staggered direction the terms (U(i,j+1) for instance) can be expressed as

Ui,j+1 = 2uN - Ui,j

We can move the 2uN to the RHS of the equation so that we do indeed get our factor 2. However the equation on the left hand side now reads (for the boundary point)
Ui,j -dt/Re(.../hx^2 +(Ui,j-1 -3Ui,j)/hy^2)

Is this indeed the correct formulation?

Quote:
 Now getting back to the issue with the boundary conditions for pressure. Have any of you written an iterative solver for the pressure poisson equation? Dr. Seibold's code uses a direct method to solve the the poisson equation for the pressure. I am facing issue with my iterative procedure for solving the equations. Plus, when the out flow boundary condition is partial(u)/partial(n) = 0, how can you say that P=0 or dp/dn=0 is a correct boundary condition? Cheers PP
When you have partial(u)/partial(n) = 0 this also implies partial(v)/partial(y)=0 (continuity).
starting from 2D NS
u_t +p_x = (u_{xx} +u{yy})/Re -(u^2)_x -(uv)_y
Assuming steady state and using partial(u)/partial(n)=0 and partial(v)/partial(y)=0 leave us with:
p_x = 0 or p = constant

We can choose this constant, for instance p=0 at outflow.

Can you explain what issues you are facing when dealing with the iterative solver for your poisson equation?

Regards,

Vincent

 September 6, 2011, 06:20 Pressure Bcn #29 Member   Prasanth P Join Date: May 2009 Posts: 30 Rep Power: 10 Hello Vincent, Your formulation is correct. Regarding the boundary condition issue, so you say that you are trying to specify the boundary condition for the steady state. So this cannot be used for an evolving flow or turbulent flow. u_t=0 -->steady state u_xx,(u^2)_x=0--->partial(u)/partial(x)=0 But why is u_yy and (v(u)_x)=0? I have tried both P=0 and dp/dn =0. For the former case the poisson solver takes several thousand some times even lakhs depending on the domain size. For the latter case the poisson doesnt converge and doesn't diverge either. dp/dn =0 is not a correct boundary condition when the velocity boundary condition is not dirichlet. But, I dono what is the issue with my solver. Cheers PP

September 6, 2011, 08:36
#30
New Member

Vincent
Join Date: Jul 2011
Posts: 29
Rep Power: 8
Quote:
 Originally Posted by prasanthnitt Hello Vincent, Your formulation is correct. Regarding the boundary condition issue, so you say that you are trying to specify the boundary condition for the steady state. So this cannot be used for an evolving flow or turbulent flow.
Indeed for a turbulent flow we do not expect p=0 at the outflow (this would imply a steady field). I guess an evolving flow means periodic boundary conditions, so p is not equal to zero in that case.

Quote:
 u_t=0 -->steady state u_xx,(u^2)_x=0--->partial(u)/partial(x)=0 But why is u_yy and (v(u)_x)=0?
I guess you mean to ask why u_yy and vu_y are equal to zero. In my previous post I took for granted that partial(u)/partial(n) =0.

For steady 3D flow between 2 parallel plates u=(u,0,0) we can apply Navier-Stokes.

u_x= 0
from NS: p_x = u_{yy}/Re (note v = 0)
lhs is only function of x and rhs only of y so p_x = constant

thus
p(x) = c*x+p0

We can take this p0 to be equal to zero for the outflow boundary.

Quote:
 I have tried both P=0 and dp/dn =0. For the former case the poisson solver takes several thousand some times even lakhs depending on the domain size. For the latter case the poisson doesnt converge and doesn't diverge either. dp/dn =0 is not a correct boundary condition when the velocity boundary condition is not dirichlet. But, I dono what is the issue with my solver. Cheers PP
If you prescribe \partial p / \partial n at all the boundaries then you obtain a singular matrix with a determinat of zero.
Could you explain a bit more on what steps you think your solver should make?
Could you specify all the boundary conditions and problem you are facing?

Regards,

Vincent

 September 6, 2011, 11:49 #31 Member   Prasanth P Join Date: May 2009 Posts: 30 Rep Power: 10 Hello Vincent, I am not interested in a steady state boundary condition. Evolving flow need not be periodic always. The boundary condition that you have given is for a very special case where v=0. Lets forget the issue about the boundary condition for the time being. Say, I use P=0 at the outflow. With my iterative Poisson solver it takes several thousand iteration or sometimes a few lakhs per time step. When I check the continuity equation at each of the pressure nodes its no where close to zero. I am not able to spot the bug in my code. So I would like to see some one solve the channel flow problem or any inflow/outflow problem with an iterative Poisson solver. The velocity is specified at all the other three boundaries. Since the velocity Bcn is a Dirichlet at these 3 boundaries, dp/dn=0 is the boundary condition for the Poisson(tats wat I assume). I am using fractional step method for solving the gov equations. Cheers PP

 September 6, 2011, 13:19 #32 New Member   Join Date: Jul 2011 Posts: 25 Rep Power: 8 Hi Prasanth, What did you use for the viscosity term? If explicit schem, try change to implicit one. Cici

 September 6, 2011, 14:09 #33 Member   Prasanth P Join Date: May 2009 Posts: 30 Rep Power: 10 Hello Cici, My viscous part is implicit. Even you had faced this problem while converting your Lid Driven Cavity code for a duct flow. Is your problem solved? Cheers PP

 September 6, 2011, 14:44 #34 New Member   Join Date: Jul 2011 Posts: 25 Rep Power: 8 Well, I can't say it is solved. Since I think both explicit and implicit viscosity will work for this problem. However, it turns out both converge now, but the solutions from explicit method are way off. So I can't say I fully understand what is going on.

 September 7, 2011, 00:36 #35 Member   Prasanth P Join Date: May 2009 Posts: 30 Rep Power: 10 Hello Cici, Good to hear that. Can you tell me the modifications that you made to your initial Cavity code. How are you solving your Poisson equation? Iteratively or some kind of direct Poisson solver? Cheers PP

 September 7, 2011, 13:06 #36 New Member   Join Date: Jul 2011 Posts: 25 Rep Power: 8 Hello, The modifications I made are the definition of boundary conditions and the schemes for solving the viscosity. The boundary conditions are described previously, and the scheme has been changed from fully explicit to implicit, ADI in my case. As to the Poisson equation, it is solved iteratively, just make sure that the boundary conditions for pressure are represented appropriately. Best, Cici

 September 7, 2011, 15:20 #37 Member   Prasanth P Join Date: May 2009 Posts: 30 Rep Power: 10 Hi Cici, "just make sure that the boundary conditions for pressure are represented appropriately" I did not quite get your point. Can you please explain. I am using Gauss Jacobi Method to solve my Poisson iteratively. I would like to know which technique are you using. This is just to know if the limitation is due to the iterative method that I am employing. Cheers PP

 September 8, 2011, 13:19 #38 New Member   Join Date: Jul 2011 Posts: 25 Rep Power: 8 Hi, I have employed the Gauss-Seidel methd for solving the pressure Poisson equation. For the boundary nodes, you have to take care of the coefficients of different terms according to the boundary conditions you defined. For instance, Coeff_e = 0; Coeff_w = 1/dx^2; Coeff_n = 1/dy^2; Coeff_s = 1/dy^2; Coeff_p = -3/dx^2 -2/dy^2; if the outlet boundary is defined as P = 0. Cici

 September 13, 2011, 05:22 #39 New Member   Vincent Join Date: Jul 2011 Posts: 29 Rep Power: 8 Since I have never worked with Gauss-Seidel iterative solver yet, I'll break the idea down and build it up from the start. We want to solve the poisson equation: Here I added a minus sign for convinience later on. We can rewrite this as a linear system that we want to solve: Where A is our negative second difference matrix and b is our right hand side of the first equation. In order to solve this equation iteratively, we use a preconditioner matrix P. This results in the following equation: The Gauss-Seidel idea is to use an preconditioner P=D+L, where D is the diagonal and L is the lower triangular part of A. Let's assume we want to solve the poisson equation in 2 dimensions. This implies a matrix K2D, which can be build by K2D = kron(I,K)+kron(K,I) Where K is the matrix approximating the minus second derivative. I think it is most insightful if I change the code to use a Gauss-Seidel solver for the pressure, if you have any question please don't hesitate to ask. I'll post it below since my post is too long. I do have to say it's not as fast as the direct method, but you can tune down the number of iterations (i=100). Also if you want a fast method you might want to look at multigrid (for 3D applications). As to your problem, the speed of convergence depends on the eigenvalues of the matrix M=I-P\A. For convergence they need to be smaller than 1, and values close to one will take a long time to dampen. The error e, can be calculated as: Good luck! Last edited by VincentD; September 14, 2011 at 03:00.

 September 13, 2011, 05:26 #40 New Member   Vincent Join Date: Jul 2011 Posts: 29 Rep Power: 8 Enjoy! Code:  function [Ue,P] = mit18086_GS %MIT18086_NAVIERSTOKES % Solves the incompressible Navier-Stokes equations in a % rectangular domain with prescribed velocities along the % boundary. The solution method is finite differencing on % a staggered grid with implicit diffusion and a Chorin % projection method for the pressure. % Visualization is done by a colormap-isoline plot for % pressure and normalized quiver and streamline plot for % the velocity field. % 07/2007 by Benjamin Seibold % http://www-math.mit.edu/~seibold/ % Feel free to modify for teaching and learning. % % 08/2011 by Vincent van Dijk % Changed to inflow/outflow conditions % Changed poisson solver to Gauss-Seidel %----------------------------------------------------------------------- Re = 1e0; % Reynolds number dt = 1e-3; % time step tf = 1e-0; % final time lx = 1; % width of box ly = 1; % height of box nx = 30; % number of x-gridpoints ny = 30; % number of y-gridpoints nsteps = 10; % number of steps with graphic output %----------------------------------------------------------------------- nt = ceil(tf/dt); dt = tf/nt; x = linspace(0,lx,nx+1); hx = lx/nx; y = linspace(0,ly,ny+1); hy = ly/ny; [X,Y] = meshgrid(y,x); %----------------------------------------------------------------------- % initial conditions U = zeros(nx-1,ny); V = zeros(nx,ny-1); % boundary conditions uN = x*0; vN = avg(x)*0; uS = x*0; vS = avg(x)*0; uW = avg(y)*0+1; vW = y*0; uE = avg(y)*0; vE = y*0; %----------------------------------------------------------------------- Ubc = dt/Re*([2*uS(2:end-1)' zeros(nx-1,ny-2) 2*uN(2:end-1)']/hy^2+... [uW;zeros(nx-3,ny);uE]/hx^2); Vbc = dt/Re*([vS' zeros(nx,ny-3) vN']/hy^2+... [2*vW(2:end-1);zeros(nx-2,ny-1);2*vE(2:end-1)]/hx^2); fprintf('initialization') Lp = kron(speye(ny),K1(nx,hx,1,3))+kron(K1(ny,hy,1,1),speye(nx)); perp = symamd(Lp); Rp = chol(Lp(perp,perp)); Rpt = Rp'; % Changed to use Gauss-Seidel for the poisson equation D = diag(diag(Lp)); L = tril(Lp)-D; Pr = D+L; % The preconditioner matrix M = speye(nx*ny)-Pr\Lp; % The multiplication matrix Lu = speye((nx-1)*ny)+dt/Re*(kron(speye(ny),K1(nx-1,hx,2,1))+... kron(K1(ny,hy,3,3),speye(nx-1))); peru = symamd(Lu); Ru = chol(Lu(peru,peru)); Rut = Ru'; Lv = speye(nx*(ny-1))+dt/Re*(kron(speye(ny-1),K1(nx,hx,3,3))+... kron(K1(ny-1,hy,2,2),speye(nx))); perv = symamd(Lv); Rv = chol(Lv(perv,perv)); Rvt = Rv'; %Lq = kron(speye(ny-1),K1(nx-1,hx,2))+kron(K1(ny-1,hy,2),speye(nx-1)); %perq = symamd(Lq); Rq = chol(Lq(perq,perq)); Rqt = Rq'; fprintf(', time loop\n--20%%--40%%--60%%--80%%-100%%\n') for k = 1:nt % treat nonlinear terms gamma = min(1.2*dt*max(max(max(abs(U)))/hx,max(max(abs(V)))/hy),1); Ue = [uW;U;uE]; Ue = [2*uS'-Ue(:,1) Ue 2*uN'-Ue(:,end)]; Ve = [vS' V vN']; Ve = [2*vW-Ve(1,:);Ve;2*vE-Ve(end,:)]; Ua = avg(Ue')'; Ud = diff(Ue')'/2; Va = avg(Ve); Vd = diff(Ve)/2; UVx = diff(Ua.*Va-gamma*abs(Ua).*Vd)/hx; UVy = diff((Ua.*Va-gamma*Ud.*abs(Va))')'/hy; Ua = avg(Ue(:,2:end-1)); Ud = diff(Ue(:,2:end-1))/2; Va = avg(Ve(2:end-1,:)')'; Vd = diff(Ve(2:end-1,:)')'/2; U2x = diff(Ua.^2-gamma*abs(Ua).*Ud)/hx; V2y = diff((Va.^2-gamma*abs(Va).*Vd)')'/hy; U = U-dt*(UVy(2:end-1,:)+U2x); V = V-dt*(UVx(:,2:end-1)+V2y); % implicit viscosity rhs = reshape(U+Ubc,[],1); u(peru) = Ru\(Rut\rhs(peru)); U = reshape(u,nx-1,ny); rhs = reshape(V+Vbc,[],1); v(perv) = Rv\(Rvt\rhs(perv)); V = reshape(v,nx,ny-1); % pressure correction rhs = reshape(diff([uW;U;uE])/hx+diff([vS' V vN']')'/hy,[],1); % p(perp) = -Rp\(Rpt\rhs(perp)); % P = reshape(p,nx,ny); c = Pr\rhs; p = zeros(nx*ny,1); for i = 1:100; p = M*p+c; end P = reshape(p,nx,ny); U = U+diff(P)/hx; V = V+diff(P')'/hy; uE = U(end,:);% + hx*Re*P(end-1,:)/dt; Will become unstable for small dt % You can remove it because P(end,:)->0 P=-P/dt; % This is the actual pressure % notice the dt missing in the pressure correction step % visualization if floor(25*k/nt)>floor(25*(k-1)/nt), fprintf('.'), end if k==1|floor(nsteps*k/nt)>floor(nsteps*(k-1)/nt) % stream function %rhs = reshape(diff(U')'/hy-diff(V)/hx,[],1); %q(perq) = Rq\(Rqt\rhs(perq)); %Q = zeros(nx+1,ny+1); %Q(2:end-1,2:end-1) = reshape(q,nx-1,ny-1); clf, contourf(avg(x),avg(y),P',20,'w-'), hold on %contour(x,y,Q',20,'k-'); Ue = [uS' avg([uW;U;uE]')' uN']; Ve = [vW;avg([vS' V vN']);vE]; Len = sqrt(Ue.^2+Ve.^2+eps); quiver(x,y,(Ue./Len)',(Ve./Len)',.4,'k-') hold off, axis equal, axis([0 lx 0 ly]) P = sort(P); caxis(P([1 end])) title(sprintf('Re = %0.1g t = %0.2g',Re,k*dt)) drawnow end end fprintf('\n') %================================================= ====================== function B = avg(A,k) if nargin<2, k = 1; end if size(A,1)==1, A = A'; end if k<2, B = conv2(A,[1;1]/2,'valid'); else B = avg(A,k-1); end if size(A,2)==1, B = B'; end function A = K1(n,h,a11,a22) % a11: Neumann=1, Dirichlet=2, Dirichlet mid=3; A = spdiags([-1 a11 0;ones(n-2,1)*[-1 2 -1];0 a22 -1],-1:1,n,n)'/h^2; happy, zhouxman and 6863523 like this.

 Thread Tools Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post h_amooie OpenFOAM Programming & Development 0 June 9, 2011 11:17 S. D. Ding Main CFD Forum 0 July 23, 2002 02:01 John C. Chien Main CFD Forum 54 April 23, 2001 08:10 ram Main CFD Forum 5 June 17, 2000 21:31 Francisco Saldarriaga Main CFD Forum 1 August 2, 1999 23:18