# Matlab code for pipe flow

 User Name Remember Me Password
 Register Blogs Members List Search Today's Posts Mark Forums Read

 October 8, 2011, 04:40 #41 New Member   Vincent Join Date: Jul 2011 Posts: 29 Rep Power: 7 I was wondering prasanthnitt, did this solution work out for you? Next it might be nice to change this code to 3D using multigrid method. Anyone up for that? Regards, Vincent

 October 8, 2011, 05:27 #42 Member   Prasanth P Join Date: May 2009 Posts: 30 Rep Power: 9 Hi Vincent, I am extremely sorry for not replying, was stuck with something and lost track of the post. I m still looking for better outflow boundary conditions. P=0 and dp/dn=-delta/dt both are working for me. Working in the sense the Poisson is converging. I have noticed one issue with the direct Poisson solver, for inflow outflow problems if we set dp/dn=0 at all the boundaries w still get some results out of it but an iterative procedure does not give a result for the same. I dono the reason for that. I m yet to implement Multigrid methods for the 2d case. But it would be a good idea to extend the code to 3d. I will be doing it in a while. I am still workin on the boundary condition issue. I need to find some for a turbulent case. People have used Outflow boundary conditions like convective boundary condition. Do any of you have any idea regarding that. Cheers Prasanth

 October 18, 2011, 12:14 I need Matlab Code #43 New Member   Mohamad Hosein Join Date: Aug 2011 Posts: 9 Rep Power: 7 hi I need matlab code that solve NS on triangular unstructured grid finite volume method. Thanks

 October 18, 2011, 14:07 #44 Member   Prasanth P Join Date: May 2009 Posts: 30 Rep Power: 9 Hi Mh. R, I am also looking for one such code. In any language its fine not necessarily in matlab. Do u have an unstructured FV code in Fortran or C or C++? If yes please share it. Cheers Prasanth

 October 19, 2011, 03:41 #45 New Member   Mohamad Hosein Join Date: Aug 2011 Posts: 9 Rep Power: 7 Hi prasanthnitt I have just fortran code that solve NS in structured grid, that written by my self, if you want it, I can send it for you.

 October 19, 2011, 05:21 #46 New Member   Vincent Join Date: Jul 2011 Posts: 29 Rep Power: 7 Why do you want to use an unstructured grid? Another option might be to use a structured grid with an immersed boundary method. Regards, Vincent

 October 19, 2011, 07:46 #47 New Member   Mohamad Hosein Join Date: Aug 2011 Posts: 9 Rep Power: 7 Hi Vincent My thesis is defined in this way. Thanks for you suggestion and attention.

 October 19, 2011, 08:46 #48 Member   Prasanth P Join Date: May 2009 Posts: 30 Rep Power: 9 Hi Mh. R, I have a structured code too, written by myself. My formulation is for an unstructured grid but for some strange reasons the unstructured results tat I get are not correct except for 2d-cavity flows. So I am working on that issue. Do let me know if you locate an unstructured incomp NS solver. Wat problem are you trying to handle with the unstructured code? Can you give the details? Cheers Prasanth

 October 19, 2011, 14:11 #49 New Member   Mohamad Hosein Join Date: Aug 2011 Posts: 9 Rep Power: 7 hi Prasanth I want solving NS with sediment transport.

 October 25, 2011, 06:26 i am new #50 New Member   David LOBON Join Date: Oct 2011 Posts: 1 Rep Power: 0 i have read the post, i am new in the topic but my question is that i want to do a dynamic state of flow inside pipe this pipe a flux external (heat transfer) and need to do 2d in the pipe and time, do you you help me ? because the scrip Benjamin is for Solves the incompressible Navier-Stokes equations in a rectangular domain...but i need inside the pipe.... i apologize but i hope to help in the foro in two week thanks

 November 12, 2011, 07:12 #51 New Member   Vincent Join Date: Jul 2011 Posts: 29 Rep Power: 7 Hi fisicas, there is no such thing as a 2D pipe, a pipe is a 3D structure. I guess 2 parallel plates is about as good as it gets as for making a 2D pipe. If your going to a 3D simulation of a pipe I recommend starting from the navier stokes equations in cylindrical coordinates. You can still treat the terms in a similar fashion as was done by Seibold. I recommend leaving the heat equation out of it until you have build and verified your 3D pipe flow code. Adding the heat equation should not be a problem however. Good luck!

 November 13, 2011, 10:46 #52 Member   Prasanth P Join Date: May 2009 Posts: 30 Rep Power: 9 Hi Vincent, Laminar Pipe flow is two dimensional. Cheers PP

 December 22, 2011, 13:17 #53 Member   Peter Join Date: Oct 2011 Posts: 52 Rep Power: 7 Hi, I have modified the same MIT18086 code but I am having some issues so I was hoping that someone in this thread who is familiar with the code might be able to help. The primary modification I did was to change the boundary condition at the lower boundary so that the velocity was equal to the upper boundary. This produces a symmetric flow with the line of symmetry at y = 0.5. Now I am trying to calculate the same flow, but only over half the domain by imposing the symmetry condition at the bottom boundary. It seems pretty strait forward but for some reason I cannot get the same answer. Does anyone here know how to do this?

 May 26, 2012, 08:29 #54 Member   Michael Moor Join Date: May 2012 Location: Ireland Posts: 30 Rep Power: 6 Hi Everyone, I am about 80% (i hope) through writing m own 2d cfd code in matlab, using a backward staggered grid, hybrid differencing, steady incompressible, and Gauss-Seidel for solving the system of equations and for convection-diffusion flow between two stationary plates... I am aving trouble however in my first iteration, with the pressure correction system of equations not converging ( it is not strictly diagonally dominant)... my question is this: Is it ok to allow the pressure corection not to converge in the first few iterations? will it get progressively more diagonally dominant as the solution develops? and hence hopefully allowing my overall solution to converge? or is this due to a poor choice of boundary conditions? NI = number of nodes in x, NJ = number of node in y i = velocity line, j = velocity line I = scalar line, J = scalar line I have attempted to use inflow outlfow bcs, the code asks the user for a stagnation pressure, which i then turn into a static one using bernoulli, and assign this to node I=1, (first scalar/pressure node), then I ask for a velocity input uin and vin, which i assign to i =2 and I=1 respectively, noting that u-cells solve from i = 3 onwards to NI-1 in the x-direction, and J=2 to J=NJ-1 in the y-direction, and v- cells solve from I=2 to I=NI-1in x, and j=3 to j=NJ-1 in y... the walls have a no-slip bc, i.e u=v=0 along all nodes in the x direction, and nodes j=2 and j=NJ-1 the outlet then, located at i=NI. says that vNI-1 = vNI, uNI=uNI-1 *(sum uin/sum uNI-1) i.e multiplying by mass out over mass in, area and density are the same. sum inflow/sum outflow. now versteeg says that the pressure gradient does NOT equal zero in the flow direction... so I have set the pressure at I= NI to zero... the pressure correction is then solved from I=2 to I= NI-1 in x, and from J=2 to J=NJ-1 in y... and all my coefficient matrices have no zero terms along the diagonal... although they are not strictly diagonally dominant... so again, if you would like to see the code i would love to put it up somewhere ( where I'm not sure as it has function files which might make it complicated to read if i put it in piece by piece here), and my questions are finally: 1. Are my boundary conditions correct? I have cut links in the momentum equations to boundaries and included the flux term from the inlet into the source term on the rhs of [A][x]=[b], and allowed uw=uw* in the pressure correction equation at the inlet, and then made no other adjustments to the pressure correction equations 2. If what I have done so far is correct, is it acceptable to use the unconverged values of the pressure correction in the solution and allow it to run it's course, or does the TDMA ( Tri - Diagonal Matrix Algorithm) solve for not strictly diagonally dominant matrices? Apologies for all the info, I hope that it is clear and you smart chaps can help me out!! Best Regards, Michael 2.

 October 2, 2012, 09:16 Special treatment of pressure correction based on continuity conservation in a pressu #55 Member   Join Date: Mar 2009 Location: Istanbul, Turkiye Posts: 45 Rep Power: 10 Hi everyone Please look this paper for the implementation of pressure correction equation boundary conditions. It is very important for algorithms you are trying to develop. SPECIAL TREATMENT OF PRESSURE CORRECTION BASED ON CONTINUITY CONSERVATION IN A PRESSURE-BASED ALGORITHM http://dx.doi.org/10.1080/10407790190053842 Nice work to you

October 3, 2012, 06:11
#56
Senior Member

Join Date: Aug 2011
Posts: 251
Rep Power: 8
Hi guys,
I just followed your discussion and here are my comments.
I guess as you pointed out that the main reason your codes do not work is due to the way you implement your BC and also what kind of BC you use especially for the pressure.

Quote:
 Originally Posted by prasanthnitt 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

dp/dn=0 is valid even if you do not use Dirichlet BC. Just consider that you use the same BC for U* (predictor step) than for U(n+1) corrector step. Then it leads also to dp/dn=0

For incompressible flows based on fractional step method (projection method, Chorin method) or SIMPLE like methods. dp/dn=0 is the right BC you should use for your poisson solver. Sometime extrapolation of pressure on the BC gives the same thing.
I know that some people claim that P=0 should be the right BC but for incompressible flows with such pressure-velocity coupling algorithm I'm a bit doubtfull. Anyway dp/dn=0 does work !
However as everyone know and as Vincent pointed out, a Poisson solver with such BC on all boundaries leads to a singular matrix.
Then you should use some specific recipes which consist in specifying one pressure node...a lot of posts on this site are concerned by this aspect..It is also not straightforward how to procede...
However some linear system solvers work well even with dp/dn=0 on all boundaries (sip solver, Jacobi, Gauss-Seidel, SOR)
For the later you just have to loop this way:

do while (res> epsilon)

loop on J
P(0,J)=P(1,J) here is included dp/dn=0
P(NI,J)=P(NI-1,J)
loop on I
P(I,0)=P(I,1)
P(I,NJ)=P(I,NJ-1)

do J=1,NJ-1
do I=1,NI-1

P(i,j)=......

end do
end do

end do

Some solvers do not admit such BC (dp/dn=0) and then they fail at least they do not converge.
Then change your solver. You can either do some stuffs in the Krylov subspace but it's a bit complicate and beyond my skills.

If you do not use FV method (but sometime it is also needed) you should correct your outlet velocity in order that the mass flow rate between inlet and outled is identical. It's a very common procedure which has been already mentioned in this discussion...

So to conclude, use dp/dn=0 with SOR or Gauss Seidel as indicated below, with velocity correction based on mass conservation and it should work.

 April 3, 2015, 06:57 FVM using SIMPLE (flow through parallel plates) #57 New Member   deewakar Join Date: Dec 2011 Posts: 12 Rep Power: 7 Hello all, Is this thread still active. I need help regarding finite volume method based code on Staggered Grid (SIMPLE method) for flow through parallel plates. I need some guidance regarding boundary conditions. Looking for reply thanks

 April 3, 2015, 07:01 #58 Member   Prasanth P Join Date: May 2009 Posts: 30 Rep Power: 9 Post your query

April 3, 2015, 08:53
here is my query
#59
New Member

deewakar
Join Date: Dec 2011
Posts: 12
Rep Power: 7
Hi

I am facing the problem with boundary conditions. I will briefly describe my problem formulation followed by query.

I have a 2D domain. nx and ny are number of nodes along x and y direction. The approach is FVM and grid is staggered backwards.

U (x velocity nodes) is defined at nodes 1,nx along x and 1,ny+1 along y.
V (y velocity nodes) is defined at nodes 1,nx+1 along x and 1,ny along y.

The pressure nodes are nx+1 and ny+1 along x and y respectively.

(Please find attached snap of the same).
The equations being solved are
Quote:
 a(J,i)*u(J,i)= a(nb)u(nb) -(P(J,I)-P(J,I-1)) .
% here nb means neighboring nodes.
% here I have written a(J,i) instead of a(i,J) to be consistent with matrix notation that i is along x , hence column number and J as row number.

here a(i,J) are the coefficients calculated based on type of scheme in terms of
Quote:
 (F= rho*u*area and D= (mu*area)/(x(J,i+1)-x(J,i)) ).
similarly I have the equation for v.

I start with initial guess u, v, p and solve the above equations to calculate u*,v*.

Here the coefficients are calculated using u and v assumed initially.
once I have these , I solve the pressure error equation

Quote:
 ap(I,J)* perror(I,J) = a(nb)*pprime(nb) + mass deficit term.
and then correct the pressure, velocity using perror terms.
The boundary conditions:
for velocity : u(:,1)=uinlet;

Quote:
 at outlet du/dx=0, u( :,nx)=u( :,nx-1) at walls u_wall=0, u(1,: )=-u(2,: ) and u(ny+1,: )=-u(ny,: )
for v velocity, v(1,: ) and v(ny,: )=0 due to wall, vin=0, which means
Quote:
 v( :,1)=-v( :,2)
v, similarly at v
Quote:
 at outlet =0, v(:, nx+1)=-v(:,nx )
P boundary conditions.
Quote:
 at exit P = 0. P(:,nx+1) dp/dx = at inlet . Hence P (:,1)=P(:,2)
In the perror equation, I am solving error for Pressure nodes from 2:nx. No error at P(:,nx+1) .

Since u velocity is known at the inlet, I am assuming error to zero at P(:,1).
The problem is after few iterations it diverges like any thing. gives NaN.
I have also attached picture of the equations that are being actually solved.
Attached Images
 SImprimante15032514521_0001.jpg (55.3 KB, 17 views) 14.jpeg (27.0 KB, 8 views)

 April 3, 2015, 09:25 #60 Member   HouKen Join Date: Jul 2011 Posts: 66 Rep Power: 7 hi everyone, I am houkensjtu, who've also been working on this problem. Its been a while(well, 3 years...), I just was too busy to put an eye on this and I hope I could help you guys all(hopefully its not too late.) Also my original matlab code is lost, but I will try to somehow reproduce the program, in matlab, python or C(I dont know yet). So for those who still working on the problem, please send me your description of your problem, to my email: houkensjtu@gmail.com happy hacking best regards,

 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

All times are GMT -4. The time now is 13:35.

 Contact Us - CFD Online - Top