Correct way to fix reverse flow at outlet
 

Hello everyone,

I'm running a low Re LES case, you can see the geometry in the attached image.

After 4s I started to have backflow from the outlet patch into my domain which made my case unstable. I was using zeroGradient for U and fixedValue 0 for p at the outlet.

I have been searching for a solution and I have found the following:

1. Change boundary conditions at outlet. Here I found two options:
   • inletOutlet for U and fixedValue 0 for p.
   • pressureInletOutletVelocity for U and totalPressure 0 for p.
   
   If there is reverse flow both options change behaviour from zeroGradient for U to not allow flow to enter the domain from the outlet. I have generally read that the second option is preferred although I don't understand the reason why.
2. Modify the solver (pimpleFoam) to have a constant K multiplying the viscosity in the whole domain. Then use setFields to define K=1 in the whole domain except for a layer of cells (~5-10) at the outlet where K=100. This way no reverse flow will be carried upstream.

So my question is, which is the correct approach? Do any of you have experience with either one to solve the issue of reverse flow from the outlet into the domain?

Thank you!

This is a typical problem which occurs with zeroGradient at output - if there is no zero Gradient in reality. The output is much to close to the input and the turbulent area to establish this. If you force the simulation to have it, you get false backflow.

This problem is not easy to resolve. You may consider setting a fixed velocity, and hoping that the symmetry works and the flow through both arms i equal. This fail if there are instability / vorticity near the output. It may also give problems with conservation of mass.

Or you ensure that the flwo at the output HAS nearly zeroGradient in U. The easiest way is to extend the arms to a length which is sufficient for level out U.

So in your view, none of the solutions I proposed is valid. Could you elaborate on why?

The problem with extending the distance to the outlet is that the number of mesh elements would increase and make my simulation run for even longer. I suppose I can try if there is no other viable option, but I'd rather not.

Hi, I am wondering if you got a solution to your problem as I am encountering the same one. Thanks.

Hi Bodo,

I solved the porblem without extending the outlet which would mean more cells and more computational time.

I would try a couple of things, in this order:

1) Use pressureInletOutletVelocity for U and totalPressure 0 for p in the outlet

2) If you don't care about the results in the outlet but you would like stability of your simulation use the outlet buffer zone region with an added constant multiplying the turbulent term in the momentum equation. For this you will need to modify the U equation and createFields file and recompile your solver. Then I recommend you to use funkysetfields in swak4foam to establish this constant as 1 in the whole domain and linearly or quadratically vary from 1 to say 100 in a region near the outlet. If you have little experience with OpenFOAM and not know how to do this, I could upload the solver.

3) If either of them (or doing both at the same time) fail, try with an advective outflow condition

If all of this fails, which I believe is unlikely you will need to extend the outlet.

Hope this helps!

Hi,

Personally, I would use the fvOption to introduce some artificial acceleration at the last outflow boundary cell to ensure all the characteristics are going outward.

Thanks,
Rdf

Hi Randolph,

Thanks for your input, I've never hear of that method. Can you elaborate? Maybe share an example of how to set up the fvOptions?

Hi,

This is not a new method. It is quite popular in DNS/LES community. I remember there is an fvOption in the newest version of OpenFOAM to apply local acceration. You can just add a source in your momentum equation at the outlet boundary cells.

Thanks,
Rdf

Dear Dreoasteh,

Thanks a lot for your reply. Unfortunately, I have a little experience with OpenFOAM, and I would greatly appreciate it if you can show me how to do the second point.

For your reference, I am using twoLiquidMixingFoam solver to simulate the mixing between two fluids. Please find attached document for some details on the problem (i.e. schematic and few screen shots of the outlet results).

Thanks for your time and I Looking forward to hearing from you.

Best Regards

Dear randolph,

Thank you. Would you please refer us to a place/website etc. where this method was applied so that we can try to mimic the methodology to our solvers? If you do not know any, would you please write some steps on how would we edit the fvOptions. I am using OpenFOAM 6.

I apologize if I am asking for a lot. Thanks again.
I look forward to hearing from you.

Best Regards

Hi Bodo,

Sure, follow this drive link to the modified solver. As you can see I only added a new field (nuk) in createFields.H and then in UEqn.H multiplied the last term of the equation by this new field nuk.

This field will not be convected or diffused. It stays the way you set it up from the time at which your simulations starts. If it has a value of 1 it does not change the simulation. If it has a value greater than 1 it increases viscosity (which is what we want to avoid back-flow which leads to instability and case blow-up). You will need to set it up as any variable, in the 0 folder.

I don't know if you compiled your own solvers before so I'll give you a more detailed explanation.

So what you want to do is download the file, unzip, move it to wherever you like (run folder is a typical choice) and with openFOAM-6 environment active go inside the twoLiquidMixingNukFoam folder and execute to compile.


This will create a solver application inside $FOAM_USER_APPBIN which you can execute as any other solver.

This will be equivalent to the original solver but, as mentioned you need to include a new file in the 0 folder to represent this new variable. A simple example (which will no impact your simulation) is also uploaded. (Note: I had to add a .txt format because this forum is old, but just remove the extension).

Then you can either use setFieldsDict to select the region you want and increase the value of nuk (generally 100 at the outler works fine) or use swak4foam.

To install swak4foam for OpenFOAM-6 follow Cyrille Bonamy's post reply. This should work fine but ask if you have any issues. This library enables you to apply mathematical expressions to openfoam utilities (i.e. set nuk varying linearly from 1 to 100 in a region defined by a box).

The utility you want to use will be called funkySetFields. You will need a funkySetFieldsDict, which I uploaded. This is the exact file I used for my simulation for the geometry seen in my first post. (Again, remove the .txt extension after downloading).

I think the usage is understandable but I will walk you through what each thing does, just in case.

create_nuk will create the nuk variable in the time directory (just in case it is not present or has weird value). It will initialize it will 1 in the whole domain and no dimensions because it is not a physical variable. This ensures that nuk is not affecting you simulation outside your specified domain.

set_nuk_top will modify the file nuk that create_nuk just made and it will replace the previous value of 1 by an expression that makes nuk vary linear depending on the value of the x coordinate (nuk=m*x+c). At the start of a the region, nuk=1 and at the end of the region (closer to the outlet) nuk=100. With these two pairs of values (x0, nuk(x0) and x1, nuk(x1)) you can make your own linear expression. Then use "condition" to define the region for which "expression" will be applied. In my case this was the x smaller than the start of my region and y bigger than the value that delimits the "upper outlet".

set_nuk_bottom does the same but for the "lower outlet". If I would have not divided it into two parts for my geometry and made only one single set_nuk instead, some parts of the inlet (which is in between the "two outlets") will have nuk>1 which I did not want - obviously.

After setting the funkySetFieldsDict to you liking just execute:



which will modify your nuk to avoid back-flow from your outlet which is what is making your case blow up.


I strongly recommend the swak4foam approach because it is a great tool to have in your openFOAM arsenal although it takes a bit longer to learn.

On a side note, I feel that you should add some more cells in the boundary layer, it seems under-resolved.


I hope this explanations helps you set up and solve your issues. Is you have any question feel free to ask, I'll be around.

Hello Dreoasteh,



Thank you for the detailed response. I am working on it.

However, I have a quick question, how would this methodology be applicable if the flow is laminar?

Thanks again and I look forward to hearing from you.

Best Regards

Hi Bodo,

I believe the method will be the same. You can apply this method whenever you encounter backflow and you can afford to lose physical accuracy in the outlet (i.e. you just want the outlet to extract the flow from the domain but you are not taking "measurements" there).

The only issue is that in a turbulent flow the appearance of eddies makes it more likely that you will encounter backflow than with a laminar flow. This implies that maybe your problem will not be solved by this method. Anyhow the test is simple: follow my comment and see what happens. At the very least you will learn about new things in the process.

Hi Dreoasteh,

Thanks for the reply.

I just realized something regarding my simulations.

For the alpha boundary conditions, I use zeroGradient or inletOutlet. Initially, it works fine since only single phase exists at the outlet boundary. However, when a mixed phase reaches the outlet boundary, I get some vortices (see attachment).
In the zero directory, I have three files only; p_rgh, U and alpha.phase1.

It seems that the vortices I get are due to inconsistent boundaries for alpha since it is only set for a single phase - alpha.phase1 (not the mixed phase).
I would appreciate your opinion about the above.

For those using inteerFoam or interMixingFoam: it is a common issue as soon as the phase with higher density touches your atmospheric boundary condition, as the pressure p_rgh is much higher within the phase of higher density and pressureInletOutletVelocity will increase U incredible here to reduce pressure to meet the totalPressure condition. In debrisInterMixing and interTempMixingFoam (both based on interMixingFoam) I apply turbulence->divDevRhoReff(rho*nuk, U) in UEqn and at the end of alphaEqns.H I do:



alpha.mixingPhase1 /= nuk;


alpha.mixingPhase2 /= nuk;


alpha.nonMixingPhase = 1.0 - alpha.mixingPhase1 - alpha.mixingPhase2;


(I do not care much about the non-mixing phase which is typically air in my mountain torrent applications, so I optimized alphaEqns.H in general to a good tradeoff between conservation of mass and bounding of concentrations to 0 < alpha < 1)