Implementation of Periodic Boundary Condition
 

Hi all,

I am going to implement periodic boundary condition in my own code. FYI, I am modeling a 2D duct with two wall mounted objects at top and bottom walls with a code that works with both SIMPLE and PISO algorithms. Does anyone have an idea how to get started? I've already done that in FLUENT but I don't know how to implement it in my code. Suggestion of papers or CFD books in this area are highly appreciated.

Thanks

Well, just connect your cells / grid points in a periodic way, i.e. in your grid management logic, and you are done!

Thanks cfdnewbie, but could you explain more or refer me to a detailed source?

I'm sorry, I don't have a good source for that. Just imagine how you would do it in one dimension with a traveling wave, that should show you how it is done...what do you find confusing about it?

I just don't know how to start!

well let's say your grid (in 1D) goes from 1 to N, which denotes the number of cells or grid points. So the periodic neighbor of cell 1 is N, and vice versa. (the other one being 2 and N-1, of course). So just when you compute your fluxes or your derivatives from your stencil, you use the periodic neighbor at the boundaries.

Did that clear it up a bit?

As stated before it is a matter of connecting your inlets and outlets as neighbours in your code. However, if you want to specify a certain mass flow you need to be able to determine the appropriate pressure drop to drive such a flow. Similarly for a constant wall temperature; you need the temperature change across the domain. A good reference for that is

"Use of Streamwise Periodic Boundary Conditions for Problems in Heat and Mass Transfer", J. Heat Transfer, Volume 129, Issue 4, 601.

I use a similar methodology in my 3D unstructured code.

It makes sense. At least I think that I have something to start with. I will be back to you in case I face any problem.

Let's recall the notation of cfdnewbie.

The way you can connect your boudary cells in a periodic way is very simple.
You just have to specify for every variables involved in your simulation just set:
PHI(1)= PHY(N-1)
PHI(N)=PHI(2)

However it's true as cdegroot mentioned it that if you have a infllow mass rate, you have to specify as source term of your Navier-Stokes equation a pressure drop which will verify the mass rate you want to impose...

I am pretty confused. As long as I am model a square (2D cube) in a duct flow, I have no problem and the solution converges. For this case I have velocity inlet and outflow boundary conditions. (outlet velocity is corrected in any iteration to satisfy continuity). But when I start periodicity as you mentioned, it gets diverged rapidly. I set those two above mentioned equalities for U, V and P; however I am dubious about what happens to my outflow after setting this condition! Is there probably any trick that I am missing?!

It sounds like your problem might be that the physics are not periodic. Do you want to model it periodic in the direction of the flow? As CDE stated, there is a pressure drop that is driving the flow. If all your variables are not periodic, then you need to do adjust them to be periodic. After all, you want to conserve mass, momentum and energy, right?

The single cube problem is not periodic but I want to model a chain of cubes in a line so that any cube (or block) is in the flow downstream of the previous block (and not in fully developed flow). Thus, the problem is periodic in the direction of the flow. I set periodicity right before solving velocity and pressure. Am I doing it correct?

I'm not sure I understand how it is periodic. Probably this is my fault. What are the reasons for running this simulation with periodic boundary conditions?

I think as DeepElm stated it that if the geometry and the pattern of the flow is not periodic it may be the reason why it diverges.
I'm not convinced that it is the case in your situation.
To be precise there is periodic boundary conditions and cyclic boundary conditions (it is periodic without pressure drop).
So you have to think what is the best option in your case...

yes ! that is what I would do...

What is your velocity-pressure coupling algorithm ?
Generally Neumann boundary conditions are used for pressure equation.
What did you do for pressure in your case ? (still neumann or periodic boundary conditions for pressure too ?)

I use SIMPLE to solve velocity-pressure coupling. I use zero gradient to set pressure of boundaries.

[QUOTE=leflix;356824]...
To be precise there is periodic boundary conditions and cyclic boundary conditions (it is periodic without pressure drop).
So you have to think what is the best option in your case...

Maybe I am wrong and as leflix mentioned, I need to consider it as cyclic. The reason I consider it as periodic (or cyclic) is that the flow is repeating. I don't want to model 10 blocks because of computational cost. If I want to model the 3 middle ones (as an example), I do not have the inlet flow because it is neither uniform in velocity or pressure nor fully developed.

You definitely need to specify a pressure drop then across the domain. However, you might be seeing the simulation diverge if you aren't setting the pressure level anywhere in the domain. If you just link up the inlet and outlet and specify a pressure drop, the actual pressure level is free to go wherever it wants and certainly will diverge. What I do is choose one volume adjacent to the outlet and set the pressure to zero at its neighbour. This keeps the pressure level from going to +/- infinity.

From your underlined comment, you might probably mean to consider a reference pressure. If so, I do consider a point, but adjacent to my inlet at (2,2), and subtract its pressure correction (from SIMPLE) from pressure corrections at all nodes. Is it what you mean or not? However, I am still confused about implementing pressure drop in my code:confused:

Sounds like you are probably correctly setting a reference pressure. Monitor your residuals though and make sure this isn't the place where the simulations starts diverging. I have found though that zeroing the pressure within an open domain can sometimes cause problems. If it is closed (like a lid driven cavity problem) there seems to be no issue. I use a fully coupled method though, so I'm not sure if you'd see the same issue using SIMPLE.

Adding the pressure drop is probably the most challenging part of implementing periodic boundary conditions. It sounds like you are considering flow past a series of bluff bodies which you are considering as periodic. I think we can agree that there should be a pressure drop across this domain. Thus, it is not sufficient to consider pressure as being periodic since the pressure at the outlet will be lower than the pressure at the inlet. Below is the way I approach the pressure drop in my code.

Consider a volume "I" adjacent to the outlet and its periodic neighbour "J" adjacent to the inlet. Instead of taking the pressure at the neighbour of "I" to be P(J) take it to be P(J)+DELTAP. Similarly take the pressure at the neighbour of "J" to be P(I)-DELTAP. DELTAP can either be specified if you know the pressure drop and want to find the mass flow. If you know the mass flow, then find DELTAP iteratively using the method of Beale described in his article "Use of Streamwise Periodic Boundary Conditions for Problems in Heat and Mass Transfer".

To set the reference pressure I choose one volume adjacent to the outlet and instead of using P(J)+DELTAP as its neighbour pressure I instead use zero. Thus I am not actually forcing the pressure to zero in any of my internal control volumes, just at a "virtual" volume adjacent to the outlet. This is sufficient to keep the pressure bounded. I hope that made sense!

I'm back and feeling pretty miserable:(
There is something with your comment that doesn't make sense to me. I am working with a collocated grid (not a staggered one). In my new geometry, my boundaries are 1 as left and N as right. As you can see it in attached file (Geometry.jpg), this is a cavity like problem with one inlet at 1 and one outlet at N. At 1, the lower half is wall and at N the upper half is wall. As long as I run this geometry as a single chamber, I have no problem and you could see the result. I just solve the equations from 2 to N-1. At 1 the velocity is prescribed and at N it is derived using outflow assumption. The contours are for U (x direction) velocity. But suppose that this is periodic problem and I am going to have the same chamber (but upside down) immediately after the current chamber. Therefore, the outlet velocity profile for the first chamber would be the inlet profile for the next one (In fact, the outlet lower half U will be the same as inlet upper half U. the same story goes with V velocity but with a negative sign). As you see, my boundary condition at 1 and N are equal with each other and not with each other's neighbors. I impose the outlet velocity to inlet at every iteration. The code converges but the problem is that this new profile does not go forward and I have it and its effect only at the very first cells (You could probably see the very narrow red area in at the inlet). After that, the velocity profile is similar to single chamber (See Geometry2.jpg). I have no idea what to do!

Hmm... so this is not your standard periodic problem since the outflow and inflow are not in line with one another. I would make two suggestions:

(i) I would suggest simulating two chambers such that your inlet and outlet are in line with one another. I don't believe that the way you have it set up that you can guarantee the outflow at the bottom is equal to to inflow at the top. If you simulate two chambers the geometric periodicity will guarantee periodicity in the flow.
(ii) It sounds like you are coupling the inflow and outflow explicitly; that is you get a certain outflow profile and apply that as a Dirichlet BC at the inlet at the next iteration. Your code will converge much better if you couple them implicitly. However this much easier if you set up the problem as suggested in (i).

Thanks for your suggestions. suggestion (i) makes sense and I'm gonna apply it know. About your second suggestion, your right; I am applying Dirichlet BC which is kinda explicit. But how can I apply it implicitly?

As an example, let's say you were computing the gradient of a quantity PHI on the east face of a volume adjacent to the outflow boundary. To do this implicitly it would be (PHI(1)-PHI(N))/DELTAX, where "1" is the volume adjacent to the inlet and "N" is the volume adjacent to the outlet. So you are treating "1" like it is the east neighbour of "N" and this would be incorporated implicitly into your coefficient matrix. You would then form all of the terms in your equations in this manner, where inlet volumes are treated like neighbours of the outlet volumes and vice versa.

Can you verify that you are conserving mass? Are you solving a compressible or incompressible system? If it is compressible, then adjusting the pressure is not sufficient for convergence.

No, it is incompressible and I do check the continuity to be satisfied.

To my mind the geometry of your problem is not periodic.
If you want to impose periodic boundary condition you should consider your domain twice as it is in your simulation. Check the sketch below, I hope you will understand what I mean.... I think this is the first reason why your simulation doesn't work. Then you will have to consider the pressure drop.

That is what I was thinking as well. It seems like you definitely need two chambers for periodicity. Thanks for the sketch by the way. That was what I was trying to explain above but was too lazy to draw out.

Yes you got the point Chris ! We agree...

Thanks leflix for the attachment. But I think Patankar solved the same problem and reported his work in:
Kelkar K.M., Patankar S.V., 1987, "Numerical Prediction of Flow and Heat Transfer in a Parallel Plate Channel With Staggered Fins" ASME J. of Heat Transfer; Vol 109, pp 25:30.
cdegroot gave me the same suggestion as you did and I am going to apply it. However, just a couple of minutes ago, I finally got result from my simulation with a single channel. The change I applied was that I first let the code converge with a prescribed inlet velocity (which prescribes a specific mass flow); then I used the converged outlet velocity as the inlet boundary condition and solved the problem again from the beginning. The thing is that if I initialize my flow field with the previous solution result, the same thing happens as before, but if I initialize it as a new problem and with zero velocity and pressure everywhere, it gets converged. I do the procedure over and over until the point that my inlet and outlet velocity difference approached to zero. I also correct my outflow with prescribed flow because otherwise my velocity increases over and over. The only thing is that Re number is a little different from what is presumed.The remedy is to correct the viscosity at any iteration to retain my Re. Now I am checking my result with the result of the paper. Then I will apply two channels side by side to investigate to all the suggestions you guys offered me.

I just found out what you meant. Thanks a lot for your help.

Coefficient matrix for Finite Volume method with periodic boundary in x direction
 

Hello,

Help please!

I wrote a 2D finite volume code for a rectangular using Neumann boundary in bottom, Dirichlet boundary at top and periodic boundary in x direction. I get some results and everything seems to be fine other than by results tilted toward right boundary and I'm assuming that is because of NOT applying correct periodic boundary condition. Could you please let me know how to build a coefficient matrix which have periodic boundary in x direction?

I appreciate if anybody can help me with that.

in x direction, if we have N control volumes. start from 0, to N-1.
where 0 means the left boundary, N-1 means the right boundary.

for control volume 1 ,
if no periodic bc: we have
AP1*1=AE1*E+Sp.
the effects from the 0 boundary has been lumped into source term Sp.
but with periodic bc, we have a fictious (N-2) control volume at the left of 1 CV. we have

AP1*1=AW1*(N-2)+AE1*E+Sp.

where AP, AE , AW means the matrix coefficent .
Hope this help, if you get some problem again, please let me know.

Hi,
I have a similar problem in implementing periodic boundary condition.
I am trying to simulate a free surface flow over an inclined surface with explicit incompressible SPH and periodic boundary conditions in x-direction. The flow is driven by gravity in x direction (gsin(theta)). But I can't get the analytical solution and the velocity keeps increasing until it converges to a value which is not what I expect (as you can see in the figure). Do you have any idea what the problem could be?
I initialise the flow with hydrostatic pressure and zero velocity.
Thanks a lot