I am a newbie. The current int
 

I am a newbie. The current interFoam solver is for incompressible fluids. Is it direct to extend this solver to compressible fluids? What should I do?

Many thanks.

I have a problem with running
 

I have a problem with running interFoam. Setup is similar to the damBreake case, i.e., 2D with symmetry boundaries, and walls. Mesh is generated in ICEM then exported in STAR (3.2) format (also tried with Fluent v6). The, using mesh utilities, converted to foam. paraFoam "sees" the mesh.
setGamma was compiled and ran on the case, mesh reread done, and in paraFoam I see my mesh and vof-scalar, as it should be.
When starting the solver, after reading enviromentalProperties, I get

--> FOAM FATAL ERROR : Cannot find a cell not on a constraint boundary starting from cell 0


Function: findRefCell(const polyMesh&, const label refCelli)

in file: findRefCell/findRefCell.C at line: 79.

FOAM exiting

Is it incorrect mesh I have imported, or I have somewhere missed some reference cell values?
Thanks
Sergei

If the mesh is properly 2D, i.
 

If the mesh is properly 2D, i.e. a slab of cells 1 cell think you can use the empty patch type for the front and back rather than symmetry boundaries which will solve your problem. Look at the definition of the damBreak mesh to see how it was done there.

Dear Weller, how about my ques
 

Dear Weller, how about my question? :-) My current problem is concerning gas and heat transfer.

The VOF idea implemented in interFOam seems to be wonderful.

It is not trivial but possible
 

It is not trivial but possible to extend the interFoam solver for compressible fluids. I have derived all the equations and have a formulation that should work well but I need sponsorship before I can allocate time to do the implementation, testing, FoamX setup etc. etc.

I have a rather simple questio
 

I have a rather simple question (before going into more details).
In one of the test cases I am running with interFOAM, a two phase flow through a micro channel junction is simulated. There is a T-like junction, which means that I have one inlet boundary (water inflow with a given velocity), and two pressure boundaries (larger cavities). Is it correct to use Outlet (two in my case) boundary condition in FOAM, and to specify a constant value of pd at those outlets?

Thanks in advance

PS. under micro channel I understand channels where capillary effects become pronounced (let say, with typical dimensions of 10-100 mu)

Yes.
 

Yes.

Henry, Thanks for the answer
 

Henry,
Thanks for the answer in relation to the interFOAM boundary conditions.
I found that my problems were related to the initialization of the solver. It was not enough to define gamma=1 and a certain velocity at the channel inlet. A proper initialization of the flow (some few cell filled with "heavy" fluid, and velocity of the fluid)is required to run interFoam smoothly.
Unfortunately, I found similar advices (by Ali) on the forum after having spent some hour on starring at deltaT going down to some 1.e-10, and even lower, and finally having resolved the problems I had.

Thanks

Yes initialisation is a proble
 

Yes initialisation is a problem with the current interFoam, you cannot specify an interface (gamma 0->1 transition) at an inlet, but I have managed to find a better way of handling inlets in the interFoam series of codes for the 1.2 release which fixes this problem.

Henry, Thanks for the info ab
 

Henry,
Thanks for the info about initialization of interFOAM solver. According to my experience, there is a large number of problems where it is quite important.
May be some dummy iterations to get interface were(?) it should be. This seems however to be problem dependent.
It will be iteresting to see how you have solved the problem in FOAM 1.2.

Sergei

I have been running some test
 

I have been running some test with the interFOAM solver, and have a question related to the Courant number. How the maximum Courant number is calculated: based on all the grid cells or only in the cells adjacent to the sharp interface (i.e. were 0<gamma<1)?
Is it possible to choose between the options?

Thanks
Sergei

based on all the grid cells bu
 

based on all the grid cells but for some problems it might be useful to calculate a maximum Courant number for the region around the interface which is quite easy to do.

Henry, Thanks for the answer
 

Henry,
Thanks for the answer with respect to Courant number in interFOAM. That is actually what I suspected. At the moment, when the flow is driven by a "heavy" phase, performance is quite OK (like in dam breake similar cases, or pipe flow of a fluid "pressed" into a pipe). In tests with a few outlets with different pressure (i.e., when the air flow is also present), computation times are really huge, because time step becomes artificially small. In some cases I have got dt=1e-8, whereas based on the heavy phase it should be about 5e-5.
Henry, is there something to do about it, and do you have intention to implement this in 1.2 version. If not, may be you can give a hint how to do it (may be I will try, which is with my c++ programming skills really an issue).

Sergei

I haven't implemented anything
 

I haven't implemented anything in 1.2 that will help you with this problem but there is something you can do about it yourself. Take a look at

OpenFOAM-1.1/src/cfdTools/incompressible/CourantNo.H

If you create a copy of this file into a copy of the interFoam directory you can make changes to it and compile your own version of interFoam with those changes. You can now redefine Courant number with any kind of biasing you like, interface, phase-fraction etc. etc. But be careful it isn't only the interface-capturing part of the solution that suffers from a Courant-number limit, the PISO pressure-velocity algorithm also suffers from bad convergence if the Courant number if larger than 1. You will have to play around to find out what you can get away with.

It may be that your problem is near-steady with respect to the air-flow and you could run with a much larger Courant number and still maintain sufficient accuracy. This can be acheived using a transient version of SIMPLE rather than PISO which is more expensive but more robust and doesn't formally have a Courant number limit. I have a test version of interFoam running with transient SIMPLE but it still has some and I don't have time at the moment to finish testing it create a releasable version.

Henry, Thanks, I will look at
 

Henry,
Thanks, I will look at the sources, and experiment with it.
I have some experience with VOF in Star (piso) and Comet (simple). My belief is, indeed, that the behavior of PISO is more predictable. However, it is really time consuming (especially when dealing with flows of fluids of very different viscosity, for example). Here, SIMPLE should perform much better, may be because of a reason you have mentioned.
When simulating flows of water-in-air relevant to our problems, I would prefer PISO with control over the Courant number definition, because in those cases interface is still rather dynamic.

Thanks again
Sergei

Will it be possible to use int
 

Will it be possible to use interFoam without a second fluid, i.e. simulate one fluid with a free surface in vacuum ?
If this is possible, what should I specify as properties for the second fluid ?

Kind regards,

ROB

No that is not possible, the V
 

No that is not possible, the VOF method is a two-fluid approach not a free-surface method.

Hi Henry, I have some few que
 

Hi Henry,
I have some few questions related to Courant number in interFoam. I have looked at

OpenFOAM-1.1/src/cfdTools/incompressible/CourantNo.H

as you have suggested above. Can you give me some more hints/additional explanations about gamma scalar field definition, which I can use in

max(mesh.surfaceInterpolation::deltaCoeffs()
*mag(phi)/mesh.magSf())

Specifically, how the gamma field can be defined/called to use it in the above construction.
What exactly

interface.cGamma()

means?
Thanks in advance
Sergei

The gamma field is a volScalar
 

The gamma field is a volScalarField but you probably would rather have it interpolated onto the faces which you can do with fvc::interpolate(gamma).

interface.cGamma() simply returns the compression coefficient you specify in the fvSolution dictionary. Take a look at the interfaceProperties.[HC] files supplied with interFoam to see for yourself.

Henry, Thank you, this is wha
 

Henry,
Thank you, this is what I was looking for. With interpolation it behaves now "better".
I will try a few bias function to speed up calculations, and compare it to Star with respect to the calculation time.

Thanks again
Sergei

There are many things you can
 

There are many things you can play with in interFoam to speed it up, you could run with a higher Courant number but you may need to increase the number of gamma sub-cycles or reduce the compression effect either by reducing the cGamma or by using upwind on the compression term. You might also want to make the solver tolerances less tight but be careful that the continuity error doesn't introduce unboundedness in gamma. Depending on you case you might also find using the AMG solver on the pressure and tuning the parameters helps.

I have tried already a higher
 

I have tried already a higher Courant number (above 0.35 I am considering not wise to go), and I have used compression coefficients of 1.5 and 1.0 (may be I have to try lowers ones). I have not varied much differencing schemas (I am using Gamma schema at the moment). This is because of analogy with Star, where I found that MARS has a superior performance for VOF, and I thought that Gamma and MARS are rather similar.
I will try also upwinding.
Thanks again for your suggestions, those are very helpful.

Sergei

The performance of the code an
 

The performance of the code and of the interface capturing are quite sensitive to the choice of scheme for both the div terms in the gamma equation. While Gamma01 with a coefficient of 1 is a good choice for the normal convection term it is less clear if this is also appropriate for the compression term because NVD/TVD analysis is not valid in this case.

I am developing a new type of scheme for the compression term which does not rely on NVD/TVD to aid boundedness and will release a preliminary version with 1.2.

Hi Henry, I have been testing
 

Hi Henry,
I have been testing interFoam with different discretisation schemes, following your advices.

1. Basically, I am currently using Gamma scheme for divergence (on gamma) with a coefficient between 0.8 and 1 (probably, for div(phi,gamma) I will use 0.5). For velocity, I am using linearUpwind. This seems to speed up the computations (in my case), especially during a number of first iterations.

2. According to your advice, I have also tried AMG for the pressure. Indeed, varying parameters speeds up the computations by some 10-15 %. I have always been thinking that AMG is more effective on large grids (I have about 100k 3D test-case, which is no t that large).

3. Results seem to be similar to STAR-CD, however, computation time is still very different (factor of about 4-5 slower). The only principal difference is that in STAR I have used an option where Courant number is restricted for the interface only. (Actually, I suspect that for interface it is kept under, let say, 0.3, and for the bulk it is probably restricted by 1.).
I have tried some variants of biasing, but not with visible success. When the Courant number is not restricted in the bulk, but only at the interface, then after a while solution does not converge (it reaches 5001 iterations, which is probably limit you built into the code, and then - diverges). This is basically, what you have been warning me about.
I tried also some gradient based biasing, but probaly made some mistakes.
Now I am thinking of something simple like maximum Courant number of 0.3 for "heavy" fluid, and 0.9, for instance, for air.
Can you advise whether the following construction is correct (I would like to multiply your definition inside the max() function in CourantNo.H by)

(1.0 - Const * neg(fvc::interpolate(gamma)))

where Const is a bias constant.
The above expression is supposed to be 1, if gamma>0, and a lower value when gamma = (or <) 0.

Can I use fvc::snGrad(gamma) for calculating gradient of the vof-scalar (if it is not zero, then there is the interface)?

Thanks in advance (it is rahter long message so)
Sergei

1) I would recommend using Gam
 

1) I would recommend using Gamma with a coefficient of 1 on both terms of the gamma equation partly because in this way the scheme is TVD conformant but also because with this value it produces near-symmetric phase-fraction distributions. I have better schemes in the 1.2 release but Gamma should be finr for your purpose.

For velocity I wouldn't recommend linearUpwind, try Gamma2V 1.0. Also it would be interesting to see what the bahaviour is like if you try upwind.

What is the Reynolds number of your case?

3) is the difference between the STAR-CD and interFoam speed due to the cost per time-step or the number of time-steps required, i.e. the Courant number?

As far as case optimisation goes I can't offer much assitance unless I play with the case and get a feel for the problem. I would be happy to do this for you if you with to contact me to do so.

Hi Henry, In my case the Reyn
 

Hi Henry,
In my case the Reynolds number is not higher than 300.
It is basically multiphase flow in a microfluidic device. In a simple test it is a t-like junction of channels, which has one inlet (with a given inlet velocity for the "heavy" phase), and two pressure boundaries with different pressure values.
This means that there is simultaneously air flow from one pressure boundary to another, and water inflow from the inlet.
At the beginning typical water velocity is around 1 m/s, whereas the air velocity is some 30-35 m/s.
If I am computing time steps for typical cells, which is based on the water velocity and taking the Courant number of 0.3, it is about 1.5*10^-5 sec.
In STAR-CD I am using time step of 10^-5 sec. Typical computation lasts for 8 hours on 4 processors.
Now about the same case in FOAM. After some few iterations, the time step goes down to about 2-3*10^-7 sec.
CPU time required per time step is shorter in FOAM, but obviously not by a factor of 30.
On one CPU case is running for longer than 5 days. This would mean that FOAM solver is more than 5 times slower.
Actually, it seems to me clear that this is primerily (I hope) because of maximum global Courant number restricted by 0.3.

I have picked up one of the typical cases we have been dealing with using STAR-CD, which have also been validated expiremtally. In this way, I can see how well FOAM is performing.
Further, I am planning to look also at conjugate heat transfer and FSI (these are main areas of interest at the moment).

Thanks for your help so far. I do understand that you have quite many messages per day, and will try not to bother you much with this. If I get some progress, I will report it.

Thanks
Sergei

Look carfully at where the vel
 

Look carfully at where the velocity is causing the maximum Courant number and why. I guess there is some level of numerical instability causing spikes in the velocity field which need to be delt with. Are these spikes at boundaries? near the interface? near corners? Are they reduced or removed by using upwind? Also how good is the mesh? Try running checkMesh on it if you are unsure about the quality.

Hi, Sergei, I am wondering
 

Hi, Sergei,

I am wondering if you have made any progress?

I am working on a project involves a small 2D channel (3 mm wide channel converges to 1 mm wide channel with 45 degree convergences. This 1 mm wide channel is 2 mm long. Then, it diverges at 45 degree to 3 mm wide). Surface tension force is strong with 70 contact angle. Results of water flowing into the convergent section seemed OK, but, I got strange results when water/air interface flowing passed the divergent section.

Could you post the schemes you used? Thanks!

Pei

Hi Pei, Basically, I had with
 

Hi Pei,
Basically, I had with interFoam reasonably goed resualts - my problem was (and is) that the solver is very "slow". I have tried to modify it wit relaxed Courant number for bulk flow (<1 in bulk, <0.3 at interface), and get speed up by a factor of 1.2-2 for my cases, but it becomes somemtimes unstable (especially when bubbles are formed).

Anyway, back to your question. You problem seems to be rather standard. Can you give me some more information (boundary conditions: pressure, inflow velocity, etc.), so I can try it during weekend.
At the moment (without knowing the exact boundary conditions) it looks like you have problems with mesh.

Sergei

Hi, Sergei, Thanks for the
 

Hi, Sergei,

Thanks for the reply!

I have run the same problem with finer meshes, different div schemes - all lead to the same characteristics. Now, I am wondering whethere this is due to strong surface tension effect (smallest length scale is 1 mm with about 10 cells, that is 0.1 mm for delta X).

The liquids I am using are air and water (1000 density ratio). I am wondering how much the solution will be affected if I reduce the density raio to 100?

If you don't mind, can I email you my case file and the mp2 movies of the results? It will be highly appreciated if you can simply take a look at the schemes I used.

Pei

I have also experienced diffic
 

I have also experienced difficulties (decreasing timestep by orders) in dealing with high density ratio cases 800:1, but when I use 40:1 I have no problem. surface tension is not high, however, maybe high shear in my problem (liq jet in air), makes it difficult to solve. The schemes I'm using:

gradSchemes
{
default Gauss linear;
grad(U) Gauss linear;
grad(gamma) Gauss linear;
}

divSchemes
{
div(rho*phi,U) Gauss upwind;
div(phi,gamma) Gauss Gamma01 0.2;
div(phirb,gamma) Gauss Gamma01 1;
}

laplacianSchemes
{
default Gauss linear corrected;
}

interpolationSchemes
{
default linear;

Hi, Ali, Thanks for the pos
 

Hi, Ali,

Thanks for the post.

I have tried out the same schemes you posted using OpenFOAM-1.1. Under OpenFOAM-1.2, I used

div(rho*phi,U) Gauss limitedLinearV;
div(phi,gamma) Gauss limitedLinear01 1;
div(phirb,gamma) Gauss gammaCompression 1;

I noticed very high vecocities near the interface. I am wondering whethere this is due to spurious currents.

Is there anything in OpenFOAM VOF that attempted to reduce spurious currents? If not, there is a paper published by Y. Renardy and M. Renardy at Virginia Tech that claims that their approach eliminates spurious currents. I am wondering if it is difficult to implement it into OpenFOAM.

Pei

Hi Pei, I guess for your ca
 

Hi Pei,

I guess for your case (surface tension dominated flow), it may be good to use filter (kernel) to smooth the volume fraction a couple of times in each time step, some improvements are observed. Some kernels are reported by Doug Kothe and Co (LANL) you may want to check and implement. They are simplest things one can do to reduce spurious currents. I believe Renardy's approach is developed for geomtery based VOF methods and many other tricks to reduce these currents are often developed for a specific case or for very simple grids. The main problem is accurate determination of interface normal and curvature where Renardy's method helps. However, it may be feasible to implement in OF.

Hi Pei, Please send your case
 

Hi Pei,
Please send your case to my home address at:
s.shulepov@chello.nl
I'l try to look at it.

Basically, with other (commercial) software, which is quite similar to FOAM, we have carried out in a past very similar (to yours) simulations, and compared those with analytical solutions.
We were actually interested in interface velocity of capillary pressure driven flows, pressure distribution in both (heavy and light) fluids, etc. Results were OK.
If you are worrying about "artificial" velocities at the interface (cells) due to local curvature variations (fluctuations) - it is a different story.

Before implementing something yourself, you can try to increase then compression factor, use Gamma01 value of 1 for both cases (see advises by Henry above). I am not happy with upwind for velocities (but again, this is my experience with STAR-CD solver, and not with FOAM), albeit people of FOAM are recommending this - they know their code better.
If you want to implement your own kernel (as advised by Ali) there is a relatively simple method described in papers by Yabe and Xiao (see for instance "Description of complex and sharp interfaces...", Computers Math Applic, v.29,p15-25).
Of course references to papers by Kothe are very good, but it will cost you a bit more time, I think.
And again, if you want, I can try to modify your model and compare it with other solvers as well.

Sergei

I am not happy with upwind for
 

I am not happy with upwind for anything and do not generally recommend it except for testing and in extreme circumstances. Generally for interFoam I would recommend limitedLinearV 1 for U, limitedLinear01 1 for the normal gamma transport term and gammaCompression 1 for the compression term. If the mesh is good, the Reynolds number is low etc. etc. you should be able to use linear (central differencing) for U and if that works OK I would certainly recommend it.

Hi Henry, I know you're too
 

Hi Henry,
I know you're too busy, but could you please give a very brief hint about what gammaCompression scheme basically is and how did you derive it? thx.

Hi Pei,
The problem of OF with high density ratio is very much reduced in OF1.2. when I run the exactly same case I mentioned above with OF1.2, timestep remained constant and all went well for a density ratio of 350:1.

It's intersting to know what kind of magic Henry did to resolven the problem http://www.cfd-online.com/OpenFOAM_D...part/happy.gif. I guess it has to do with pressure-velocity coupling and also twoPhaseMixture (maybe).

If you still have problems with OF1.2, you may want to implement the kernel, but before that make sure you play with parameters in fvSolution too, you can decrease the convergence criteria and see if it improves the results or not.

Hi, Ali, Sergei, and Henry,
 

Hi, Ali, Sergei, and Henry,

Thanks for the replies and suggestions.

Sergei, I have sent you the video files and case files to your home email address.

I have used both OpenFOAM-1.1 and OpenFOAM-1.2 and tried out all the schemes (div) suggested. All lead to similar flow characterisitcs. I wish I can post these video files so that you can take a look at them. I think the results look strange. There is a possibility that this could be real, but, I do not know for sure. I have to setup some experiments to test this.

Thanks again for your feedback.

Pei

Hi, I have tested for a while
 

Hi,
I have tested for a while interFam solver of version 1.2 on different meshes and different problem setups, and have to say that VOF solver performs very well. Even on rather coarse meshes results are suprisingly (at least to me) good.
It is still much slower than say STAR, but behavior of interface is much more "physical".
Henry, it is really good code, thanks.

Further, I have a question. I am actually also ineteresting in VOF with moving meshes. I have seen some discussions on interFoam and dynamic mesh, and even tried to run some examples but without success till now (because of lack of knowledge in C++, I think).
I would like to start from a simple case: squeezing flow, i.e., 2D channel, where upper wall moves down, compressing light and heavy fluid. Wall velocity may in general be a function of, let say, time. In this way, I can look at mass, momentum, etc. balance, and there are also analytical solutions for different geometries and setups (2d, axi-symmetric).

Can somebody gibve me a suggestion on how to dynamically scale the mesh during computations in a simple way (if one exests)?

Pei, may be you have already some draft version of tutorials on moving meshes? I will really appreciate it.

Again, at the moment I am interested in a simple scaling, and not in topological changes of the domain.

Thanks in advance
Sergei

PS. If this case would work (I have tried to run examples found on discussion forum, but have problems in v.1.2), and there is an interest, I can submit a simple tutorial, as a contribution into documentation.

Hi, Sergei, I wrote a draft
 

Hi, Sergei,

I wrote a draft about 2 months ago (un-complete). I had to stop because I changed to a new job. I can email you the cases Dr, Jassak sent me and the steps to run them.

But, I do not understand why you want to scale the mesh. In the situation you mentioned above, moving the wall down (or up), you only need to add mesh layers or removing mesh layers during run time.

Pei

Hi Sergei, mesh movement is
 

Hi Sergei,

mesh movement is done through the movingFvMesh class in the movingFvMesh library ($FOAM_SRC/movingFvMesh). The currently selected mesh gets its 'move()' function called whenever the mesh has to be moved. Have a look through the ones in the movingFvMesh/lnInclude directory to find the one that behaves the most like you want it to and copy/modify it.