[chegdan]

Hello everyone,

First of all, these should be quick questions for a few of you.

1. I'm trying to import a mesh file created in gambit, exported to .msh format and imported into OF-1.4.1-dev using the fluentMeshToFoam utility. Once imported there are a few boundaries with the "_shadow" attached to the end of them. How do I handle them in OF? Right now I have a cyclic boundary called "cyclic" and a shadow bc called "cyclic_shadow". I have them both defined as cyclic base types...is this correct? Does a cyclic BC need a defined value (Uniform/nonuniform?)?

2. I would like to run a case using interFoam with periodic BC's. In order to properly handle pressure do i need to do something similar to channelOodles with creategradP and writeGradP? Will the (VOF method) interfoam solver handle this?

Thanks for your help in advance.

Dan

[chegdan]

Just thought I would respond to my own post with the solution.

1. If you are using Gambit and want periodic boundaries, just link the two sides (its in the manual), mesh the object, and then define a single bc as those two faces (3d)/edges (2d) as a type "wall". After you import the mesh then go to your constant/polyMesh/boundary file and change the boundary type of your cyclic boundary to type "cyclic". If you define your initial bc's correctly it should work.

2. The short answer to this is yes, you will need to use create and writeGradP, but I am currently having problems passing gamma to the other side of my cyclic boundary.

Just thougth this might help someone trying to use gambit and make a domain with cyclic boundaries.

Dan

[chegdan]

Originally Posted by chegdan

Just thought I would respond to my own post with the solution.

1. If you are using Gambit and want periodic boundaries, just link the two sides (its in the manual), mesh the object, and then define a single bc as those two faces (3d)/edges (2d) as a type "wall". After you import the mesh then go to your constant/polyMesh/boundary file and change the boundary type of your cyclic boundary to type "cyclic". If you define your initial bc's correctly it should work.

2. The short answer to this is yes, you will need to use create and writeGradP, but I am currently having problems passing gamma to the other side of my cyclic boundary.

Just thougth this might help someone trying to use gambit and make a domain with cyclic boundaries.

Dan

I thought i would answer my own question and another that I saw related to the topic:

I messed with this last year in OF solely because I was told I would spend my life trying to get Fluent to do this properly (thankfully for that comment my advisor has switched completely to OF since). I was able to get cyclic boundaries to work with interFoam and rasInterFoam. I remember that the problem was with the fact that gravity was lumped into the pressure term and when the liquid drop fell through the bottom and was supposed to come back through the top...the liquid drop just stopped and hung in mid air with small pieces showing at the top of the domain. it was actually pretty funny to see. The solution was to add rho*g to the rhs of the momentum equation (pd = P - rho g z. i.e. the reduced pressure is not periodic, but the absolute pressure. interFoam solves pd, not P) (in addition to all the changes to make it similar to channelOodles).
The main bit of code that was really a simple change in UEqn.H was:

   fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U)
      + fvm::div(rhoPhi, U)
      - fvm::laplacian(muf, U)
      - (fvc::grad(U) & fvc::grad(muf))
    //- fvc::div(muf*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
	==
	rho*g				//added to pull gravity out of pressure term
        + flowDirection*gradP		//added for channel flow
    );

    UEqn.relax();

of course there are a few more changes for channel flow and you need to read in the g field and value of rho, but that is the gist of it. I never changed the momentum corrector part of the UEqn.H file...but it would probably need to be changed if you use that.

Lastly, one thing that i noticed was that as I dropped a liquid drop and it traveled through the domain boundaries and appeared again on the top...it kept speeding up. I am pretty sure this was because there was no drag term in the momentum equation...something that I know for the case of the rising bubble has been studied for a while. There are numerous threads on that topic. Anyway, long answer for a simple question...yes its possible to do this in OF.

[AlmostSurelyRob]

Thanks for pointing me to this post and also for keeping it up to date. I have an implementation in Fluent for channel where I correct momentum, pressure and phase fraction. This implementation doesn't work, when I turn on gravity. I guess it's due to the reasons you outlined above, but with Fluent I have little chance of investigating the origin of the problem. Your post sounds very encouraging so I'll try to rewrite it to interFoam.

[joris.hey]

hello chegdan,
I'm arriving 2 years after the discussion but I hope somebody is still watching !

As many others, I'm trying to modify the interfoam solver to an Inter-Cyclic-OpenChannel-Foam one. I model flow over bedform when water depth scale with bedform heigth. It is happenning on steep slopes in mountains river.

Here, the gravity component in x-direction is the only energy source in the system.

My idea to begin the simulation was to start with an initial volume of water at rest, that will move forward with the action of gravity and come to a steady (water depth/velocity/pressure) when the turbulent dissipation will compensate the gravity source. (see image below at the initial condition)



I made the modification you suggested to deal with cyclic conditions for relative and absolute pressure adding the rho*g term.

I didn't made the change to keep the flow rate constant (like channelOodle) as it is not constant in my case but going from zero (initial water at rest) to the limit discharge.

However, it seems I still have a problem at the cyclic patch, water seems not to go freely back to the input... Look at the images below :

Alpha :


Pressure :


Velocity:


Where do you think it is coming from ? Should I try to implement also the constant flow rate ? I guess the bubble case is a bit the same than mine for the transcient behavior ...

Maybe I missed something when you said : "separate gravity from pd and add it back into the momentum equation" ... How can I separate the gravity from pd ? in which file ?

Also "you need to read in the g field and value of rho" is not clear for me. simply adding rho*g in Ueqn.h is not enough ?

Thanks very much !
Joris

[joris.hey]

Ok, I guess I found the solution...

In addition to add rho*g in the RHS of UEqn.h , I deleted the line:

- ghf*fvc::snGrad(rho)

in both Ueqn.h and Peqn.h, trying to get rid of the substracted gravity in p_rgh.

Is it correct ?
It seems to work at first look.

Cheers,
Jors

[chegdan]

Joris,

I haven't been in this solver for a little while and its changed a little bit from 1.4.1-dev when I did some fiddling with interfoam. I did notice that they changed the algorithm a bit and now it uses the PIMPLE method (it seemed a little bit simpler a few versions of OF ago). So, with this...I am going to just say be careful when you make quick ad-hoc changes and I cannot confirm if what you did was correct. When you do find an answer, post back that it worked for others to see.

[joris.hey]

hehe yeah maybe it is a bit violent what I did, but as I'm not really familliar with the FOAM coding...

Also I forgot to say that I changed your proposition

   fvVectorMatrix UEqn     (         fvm::ddt(rho, U)       + fvm::div(rhoPhi, U)       - fvm::laplacian(muf, U)       - (fvc::grad(U) & fvc::grad(muf))     //- fvc::div(muf*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf())) 	== 	rho*g				//added to pull gravity out of pressure term         + flowDirection*gradP		//added for channel flow     );      UEqn.relax();

into

   fvVectorMatrix UEqn     (         fvm::ddt(rho, U)       + fvm::div(rhoPhi, U)       - fvm::laplacian(muf, U)       - (fvc::grad(U) & fvc::grad(muf))     //- fvc::div(muf*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf())) 	- rho*g				//added to pull gravity out of pressure term         + flowDirection*gradP		//added for channel flow     );      UEqn.relax();

as the == sign is below in the PIMPLE loop. Don't know if it changes anything.

I definitely not converge to a finite velocity, and after 10 sec waves developed on the free surface cause the celocity (and the froud number) becomes too high...

This is maybe because I used a relatively coarse mesh and I should better use RANS or LES than laminar to model small scales dissipation...

The problem with the flow I try to model is that in reality there is a strong dissipation due to the wall friction (the wall roughness is in the same order of magnitude as the flow depth) and it seems a bit strange to model it with a "flat" curved wall. I know that there exists some wall functions in OpenFOAM, I migth have to use it...

Again, thanks for advices !

Joris

[joris.hey]

Ok, here I post the results I got with what I did above. It seems nice.

Water:
http://www.youtube.com/watch?v=wEjsM...B1BLS8joo24x_y

Velocity:
http://www.youtube.com/watch?v=Rr6aQkXUlh8&context=C3f6e459ADOEgsToPDskJj rSjp4uB1BLS8joo24x_y

Thanks for helping
joris

[13msmemusman]

hi i got errors like

UEqn.H:5:23: error: ‘muf’ was not declared in this scope
- fvm::laplacian(muf, U)
^
UEqn.H:11:11: error: ‘flowDirection’ was not declared in this scope
+ flowDirection*gradP
^
UEqn.H:11:25: error: ‘gradP’ was not declared in this scope
+ flowDirection*gradP
^
Please help he how to declare these variables

Originally Posted by chegdan

I thought i would answer my own question and another that I saw related to the topic:

I messed with this last year in OF solely because I was told I would spend my life trying to get Fluent to do this properly (thankfully for that comment my advisor has switched completely to OF since). I was able to get cyclic boundaries to work with interFoam and rasInterFoam. I remember that the problem was with the fact that gravity was lumped into the pressure term and when the liquid drop fell through the bottom and was supposed to come back through the top...the liquid drop just stopped and hung in mid air with small pieces showing at the top of the domain. it was actually pretty funny to see. The solution was to add rho*g to the rhs of the momentum equation (pd = P - rho g z. i.e. the reduced pressure is not periodic, but the absolute pressure. interFoam solves pd, not P) (in addition to all the changes to make it similar to channelOodles).
The main bit of code that was really a simple change in UEqn.H was:

   fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U)
      + fvm::div(rhoPhi, U)
      - fvm::laplacian(muf, U)
      - (fvc::grad(U) & fvc::grad(muf))
    //- fvc::div(muf*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
	==
	rho*g				//added to pull gravity out of pressure term
        + flowDirection*gradP		//added for channel flow
    );

    UEqn.relax();

of course there are a few more changes for channel flow and you need to read in the g field and value of rho, but that is the gist of it. I never changed the momentum corrector part of the UEqn.H file...but it would probably need to be changed if you use that.

Lastly, one thing that i noticed was that as I dropped a liquid drop and it traveled through the domain boundaries and appeared again on the top...it kept speeding up. I am pretty sure this was because there was no drag term in the momentum equation...something that I know for the case of the rising bubble has been studied for a while. There are numerous threads on that topic. Anyway, long answer for a simple question...yes its possible to do this in OF.

[cuikahlil]

Hi guys!

I have come 6 years late to this post, and it sure would be a long shot if you can still see this. I'm also doing a cyclic boundary interfoam where my gravity is multiplied to an angle to simulate a slope. Fluid is already in my domain as initialized from setFields, and it would collapse on to a cyclic wall and would just collapse on itself over and over again.... well at least that's what I wanted it to do. I did the suggestions in this post about separating the gravity term and removing the -gh*fvc::snGrad(rho), and was even able to compile it correctly, but it still doesn't work. OF has made a lot of changes to interFoam since the time these answers were posted so I'm not sure if the results would still be the same

My fluid simply disappears from the domain as if it passes through the bottom wall. I've checked the mesh and it seems to be fine. lowering the time step also doesn't work. If it would be of any help, my Ueqn is

surfaceScalarField muEff
(
        "muEff",
	mixture.muf()
	+ fvc::interpolate(rho*turbulence->nut())
);

fvVectorMatrix UEqn     
    (         
        fvm::ddt(rho, U)       
      + fvm::div(rhoPhi, U)       
      - fvm::laplacian(muEff, U)       
      - (fvc::grad(U) & fvc::grad(muEff))//- fvc::div(muf*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf())) 	
        ==
        rho*g				//added to pull gravity out of pressure term              
    );      
    UEqn.relax();

    fvOptions.constrain(UEqn);

    if (pimple.momentumPredictor())
    {
        solve
        (
            UEqn
         ==
            fvc::reconstruct
            (
                (
                    mixture.surfaceTensionForce()
                  - fvc::snGrad(p_rgh)
                ) * mesh.magSf()
            )
        );

        fvOptions.correct(U);
    }

   fvVectorMatrix UEqn     
(

I'd be happy to provide more if necessary.
Thank you in advance, and I really do hope you can help me with this.

[cuikahlil]

Hello again!

After a bit of code editing, the solution to my problem was provided in the codes in the final panels of this thread:

interFoam and cyclic boundary condition issue

interfoam has undergone quite a number of revisions since this solution was presented, but making it work is fairly manageable.

[hrvig]

Hello foamers,


I tried to implement the modifications into the recent interFoam solver provided with OpenFOAM6 (attached). When I simulate a simple bubble rising in a stagnant fluid in a cyclic domain I see strange behaviour in the pressure field when the bubble passes the cyclic boundary. Does anyone know more specifically how to fix this issue? I have attached the simple test case I am running.


Best regards,
Jakob

Originally Posted by cuikahlil

Hello again!

After a bit of code editing, the solution to my problem was provided in the codes in the final panels of this thread:

interFoam and cyclic boundary condition issue

interfoam has undergone quite a number of revisions since this solution was presented, but making it work is fairly manageable.

[hrvig]

The problem almost solved itself. I forgot to comment out line 35 in pEqn.H:

- ghf*fvc::snGrad(rho)

After doing so it now works flawlessly. I have attached the solver if you are interested.



Best regards,
Jakob

Originally Posted by hrvig

Hello foamers,


I tried to implement the modifications into the recent interFoam solver provided with OpenFOAM6 (attached). When I simulate a simple bubble rising in a stagnant fluid in a cyclic domain I see strange behaviour in the pressure field when the bubble passes the cyclic boundary. Does anyone know more specifically how to fix this issue? I have attached the simple test case I am running.


Best regards,
Jakob

[Bodo1993]

Originally Posted by hrvig

The problem almost solved itself. I forgot to comment out line 35 in pEqn.H:

- ghf*fvc::snGrad(rho)

After doing so it now works flawlessly. I have attached the solver if you are interested.



Best regards,
Jakob

Hi hrvig,

Thanks for sharing the solver. I am wondering if the original interFoam solver you edited was from version 6 (OpenFOAM v6).

I am having similar issue. In a simple 2D vertical channel flow, with one fluid initially filling the channel and the second fluid is pumped from below. The results are correct except for the outlet boundary, where the fluids just exit from a small potion of the outlet face (I do not use cyclic boundary condition, I have totalPressure for P_rgh, pressureInletOutletVelocity for U and zeroGradiant/inletOutlet for the phase).

I have tried to mimic the edits as you have done. It fixed the outlet issue, however, it created other issues in the computational domain.

I would greatly appreciate any assistance.

[Bodo1993]

Originally Posted by 13msmemusman

hi i got errors like

UEqn.H:5:23: error: ‘muf’ was not declared in this scope
- fvm::laplacian(muf, U)
^
UEqn.H:11:11: error: ‘flowDirection’ was not declared in this scope
+ flowDirection*gradP
^
UEqn.H:11:25: error: ‘gradP’ was not declared in this scope
+ flowDirection*gradP
^
Please help he how to declare these variables

Hi,

I got the same error message:

(‘flowDirection’ was not declared in this scope)
(‘gradp’ was not declared in this scope)

I am wondering if you figured out the way to do that properly.

Would appreciate your response.

Thanks

[Neshat]

Originally Posted by joris.hey

hehe yeah maybe it is a bit violent what I did, but as I'm not really familliar with the FOAM coding...

Also I forgot to say that I changed your proposition

   fvVectorMatrix UEqn     (         fvm::ddt(rho, U)       + fvm::div(rhoPhi, U)       - fvm::laplacian(muf, U)       - (fvc::grad(U) & fvc::grad(muf))     //- fvc::div(muf*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf())) 	== 	rho*g				//added to pull gravity out of pressure term         + flowDirection*gradP		//added for channel flow     );      UEqn.relax();

into

   fvVectorMatrix UEqn     (         fvm::ddt(rho, U)       + fvm::div(rhoPhi, U)       - fvm::laplacian(muf, U)       - (fvc::grad(U) & fvc::grad(muf))     //- fvc::div(muf*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf())) 	- rho*g				//added to pull gravity out of pressure term         + flowDirection*gradP		//added for channel flow     );      UEqn.relax();

as the == sign is below in the PIMPLE loop. Don't know if it changes anything.

I definitely not converge to a finite velocity, and after 10 sec waves developed on the free surface cause the celocity (and the froud number) becomes too high...

This is maybe because I used a relatively coarse mesh and I should better use RANS or LES than laminar to model small scales dissipation...

The problem with the flow I try to model is that in reality there is a strong dissipation due to the wall friction (the wall roughness is in the same order of magnitude as the flow depth) and it seems a bit strange to model it with a "flat" curved wall. I know that there exists some wall functions in OpenFOAM, I migth have to use it...

Again, thanks for advices !

Joris

Hi Joris,

Did you finally apply the LES model? Would you please let me know which wallFunction did you use for the bottom of your bed form?

Thanks in advance
Neshat