parasitic currents
 

Hi,

I have working with high weber number flows (capillary flows) and I am finding som important parasitic currents in the interface.

I have read that Brackbill (1992) J. Comput. Phys. 100 335-354, suggested to weight the surface tension force with the density in the cell in order to reduce this non-real currents.

I have modified this term in interFoam,

fvc::interpolate(interface.sigmaK())*fvc::snGrad(a lpha1)

multiplying by

fvc::interpolate(rho/(twoPhaseProperties.rho1()-twoPhaseProperties.r
ho2())

The parasitic currents have been appreciably reduced. I am surprised than I so simple thing has not been before implemented in interFoam. Maybe there is some drawback that I don't know?

Best regards

Hello Robert,
I agree with you, that we need a new thread for this discussion.
I have a basic question, how do you usually visualize or calculate your parasitic flows at the interface. do you do this by just performing some post processing or is there any script based on some literature to calculate this. I need to know this because I am running some lesinterfoam calculations and it is really important for me to resolve the forces properly at the interface.

I am using OF-1.5, and if the implementation that you did is specified in some literature you can also let me know the paper, I will look in to it.

I may also change the lesInterfoam solver and see if there is any major difference.

But first of all i want to know how to find out if there are any currents acting at my interface,
Hope you dont mind my basic question. I am new to the field of Multiphase flows. My Master thesis was in Incompressible flows, so i dont have much knowledge of Multiphase flows.
bye
with regards
K.Suresh kumar

Hi Robert

Thanks for sharing this information. I have been hvaing some similar issues with the parasitic currents. Please post some pictures showing the comparison.

According to you post:

Does that means that your final expression in pEqn.H for phi now lokks like this:

phi = phiU +
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(g amma)*
fvc::interpolate(
rho / (twoPhaseProperties.rho1()
-twoPhaseProperties.rho2()
)
- ghf*fvc::snGrad(rho)
)*rUAf*mesh.magSf();

Thanks once again for sharing the information.

Best Regards
jaswi

Hi,

kumar, the paper by Brackbill is, as far as I know, the first on parasitic currents. You can read also Harvie et al. Applied Mathematical Modelling, 30 (2006) 1056-1066.

jaswi, not exactly, since I am using OF 1.6 and ghf is no more used. The modified pEqn.H is as follows:

phi = phiU +
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(a lpha1)*
fvc::interpolate(rho/(twoPhaseProperties.rho1()-twoPhaseProperties.r
ho2()))*mesh.magSf()
+ fvc::interpolate(rho)*(g & mesh.Sf())
)*rUAf;

Attached you can see the velocity vectors of one experiment, made with paraview. In the right side, with density correction, the parasitic currents have been appreciably reduced.

Hope it can help you.

with regards

Robert

Dear Robert

Thanks for the answer and the pictures.
They do show a substantial decrease in the spurious currents.
I just can't wait to try that for my case :-)

BR
jaswi

Thank you, Robert. This one was a really important notice (why I did not found it out by myself?). I am also dealing with low Capillary number flows and my OF model is suffering from high unphysical velocities close to interface. I will try modified model and see what would be the difference in my case.

One drawback here is that now it will take more time to compute. How much? That is an another question.

By the way, am I right that you did changes both in pEqn.H and UEqn.H ?

BR
Denis

Hi, Denis,

actually I have only modified the pEqn.H, since the Ueqn.H is only calculated if momentumPredicor is "true", which is not normally the case.

Hope it will be useful for you.

Best regards

Robert

I got your point. But to be on a safe side I modified both.


Ok, I got some intermediate results. It seems that such a simple modification completely cured that annoying effect when there is unwanted high velocities near the interface and now my OF model is much better follow the experimental data.
Thank you, Robert one more time for posting this idea here.

I do not know how, but I think this idea should be somehow more widely known, because, I think, a lot of people deal with capillary flows and might face the similar problem.

has anybody applied this method for static problems?
 

today, I made similar changes to interFoam, but it seems that it produces higher parasitic currents in my case which is kind of calculating capillary pressure in a bubble trapped in a closed micro-scale 2D shape with some convex edges (all boundaries are no-flow ->static problem). it seems that the term "fvc::interpolate(rho/(twoPhaseProperties.rho1()-twoPhaseProperties.rho2())" is greater than one and it causes higher surface tension force, and higher parasitic currents in my static problem.
has anybody applied this method for similar static problems?

I have mistyped the term. The suggestion of Brackbill is rho/<rho>, i.e.,

fvc::interpolate(2*rho/(twoPhaseProperties.rho1()+twoPhaseProperties.rho2 ())

This term has to be 1 in the interphase, and it has to reduce the force in the lighter phase. It can be greater than one in the heavier phase, but it should not increase the parasitic currents.

Robert

Thanks for the corrections.
It is noticeable that the previous version also worked fine for the gas-liquid flow as it has value of 1 in the heavier phase and at the interface and a very small value in the lighter phase.

Hi all,

I made the same modifications to a 2D static bubble problem with the result, that an original round bubble gets octagonal!! So the results get worse! Right now we performing some systematic test.
The parasitic current problem is not that simple, since they depend on the physical problem (density ratio, viscosity ratio, We, Re, ...) and the numerical implementation of the VoF method.

Regards Jens

Hi Jens,

since one year passed, do you have any results of you investigation to share?

Best Regards,
Ilya

Hi all,
i also made the same modification in pEqn.H, but i still get very high spurious velocity at the interface. if you remember, since one year passed, coul you please post your fvScheme and fvSolution. and what were the properties of your fluids?

to have a look at my problem you can read here: http://www.cfd-online.com/Forums/ope...tml#post308969

Thanks

andrea

Hi,

thanks for posting this. I implemented the change in OpenFOAM 2.1.0, but it seems to have no effect here.

As a benchmark, I am using a sessile droplet (d = 1mm, semi-sphere) on a flat plate. After settling down due to the distortions caused by the initialisation of the alpha1 phase on the rectangular mesh (about 0.03s), the drop starts to move around over the plate. I suspect this behaviour is due to spurious currents, which I hoped to prevent with the measures described above.

I will try to verify this result once more.

We have solved the problem of spurious currents in our modified version of interFoam code, you can find the details in a paper (two-phase flow at low capillary numbers / small scales) being published in JCP:

http://dx.doi.org/10.1016/j.jcp.2012.04.011


The Sharp Surafce Force part, described in the paper, is a one line code which any one can easily implement. Smoothing and filtering are also proposed in the paper.

The paper is written for Cartesian meshes, but most parts are general, and now we are using the code for unstructured meshes with very few modifications / change of discretization algorithms. Hopefully these modifications will be published soon as well.

Hello,
Can you send me the paper to pablodecastillo at gmail dot com or give here a few hints??

Thanks

As a quick reply, the Sharp(er) Surface Force(SSF) formulation is as follows:

in interFoam/pEqn.H change
fvc::snGrad(alpha1)
to
fvc::snGrad( (1.0/(1.0-2.0*C)) * min( max(alpha1,C), (1.0-C) ) )

For Static cases use C=0.47-0.49, for dynamic simulations (moving interfaces) use C=0.15-0.20, roughly speaking. (the definition of C here is different from the paper)

Filtering is as important as sharpening the capillary force but it has a long story. I will be back to you later by email.

Thanks Ali,

I will try tomorrow morning, and i will coment to you how it goes.

You linked to your email instead of doi.

I found it here:

http://www.sciencedirect.com/science...21999112001830

sorry, apparently the forum don't like the copy and paste approach !

fixed

I implemented an open source version of the solver based on Ali's paper: http://code.google.com/p/interfoamssf/

It's still buggy (the time step has to be constant, and the dambreak results are not good yet), but I decided to release it to speed up development. Contributors are very welcome.

- Anton

Hi, Anton,

I am wondering if you are willing to share your revised solver so that I can test it out?

Thanks!

Pei-Ying
phsieh2005@yahoo.com

Hi,

just write

hg clone https://code.google.com/p/interfoamssf/

in a shell and you can download it.

best


andrea

Thanks a lot Andrea!

Pei-Ying

hi,

I am wondering if anyone has successfully compiled the interfoamssf code. I got the following error:

SOURCE=interFoam.C ; wmakeScheduler g++ -m64 -Dlinux64 -DWM_DP -Wall -Wextra -Wno-unused-parameter -Wold-style-cast -O3 -DNoRepository -ftemplate-depth-100 -I/home/phsieh/OpenFOAM/OpenFOAM-2.1.x/src/transportModels/twoPhaseInterfaceProperties/alphaContactAngle/alphaContactAngle -I/home/phsieh/OpenFOAM/OpenFOAM-2.1.x/src/transportModels/incompressible/lnInclude -I/home/phsieh/OpenFOAM/OpenFOAM-2.1.x/src/finiteVolume/lnInclude -I/home/phsieh/OpenFOAM/OpenFOAM-2.1.x/src/transportModels -I./interfaceProperties/lnInclude -IlnInclude -I. -I/home/phsieh/OpenFOAM/OpenFOAM-2.1.x/src/OpenFOAM/lnInclude -I/home/phsieh/OpenFOAM/OpenFOAM-2.1.x/src/OSspecific/POSIX/lnInclude -fPIC -c $SOURCE -o Make/linux64Gcc46DPOpt/interFoam.o
interFoam.C: In function ‘int main(int, char**)’:
interFoam.C:93:21: error: ‘class Foam::pimpleControl’ has no member named ‘start’
interFoam.C:93:62: error: no ‘operator++(int)’ declared for postfix ‘++’ [-fpermissive]
interFoam.C:152:42: error: ‘class Foam::pimpleControl’ has no member named ‘nCorr’
In file included from interFoam.C:154:0:
pEqn.H:41:69: error: no matching function for call to ‘Foam::pimpleControl::finalInnerIter(int&, int&)’
pEqn.H:41:69: note: candidate is:
/home/phsieh/OpenFOAM/OpenFOAM-2.1.x/src/finiteVolume/lnInclude/pimpleControlI.H:80:13: note: bool Foam::pimpleControl::finalInnerIter() const
/home/phsieh/OpenFOAM/OpenFOAM-2.1.x/src/finiteVolume/lnInclude/pimpleControlI.H:80:13: note: candidate expects 0 arguments, 2 provided
/home/phsieh/OpenFOAM/OpenFOAM-2.1.x/src/finiteVolume/lnInclude/readTimeControls.H:38:8: warning: unused variable ‘maxDeltaT’ [-Wunused-variable]
make: *** [Make/linux64Gcc46DPOpt/interFoam.o] Error 1
phsieh@hillsdale:~/OpenFOAM/phsieh-2.1.x/solvers/interfoamssf
---------------------------------------------------

I compiled this under OpenFOAM-2.1.x (latest git pull) on OpenSUSE 12.1 64 bit.

Thanks!

Pei-Ying

I'm using OF-2.0.x. The PIMPLE controls have changed a bit in version 2.1.x.

HI, Anton,

Thanks for the reply!

Under OpenFOAM-2.1.x, I copied pimpleControl.C pimpleControl.H, and pimpleControlII.H (these 3 codes were copied from OpenFOAM-2.0.x source) into interfoamssf folder, and the solver compiled without any error.

However, when I ran the damBreak test case, I got several errors. I then, added
-----------------------
"alpha."
{
solver smoothSolver;
smoother GaussSeidel;
tolerance 1e-11;
relTol 0;
nSweeps 1;

maxIter 100;
nLimiterIter 20;
maxUnboundedness 1;
CoCoeff 0.5;
}
----------------------
into fvSolution, and
default Gauss linear; under divSchemes in fvSchemes.

The case then ran upto time - 0.4. But, when I plotted alpha1 using paraFoam, the liquid/air interface did not even move. So, there is something wrong. Any suggestion?

Once I can get the dambreak test case to run correctly, I am planning to run something with more complex geometry and not so nice-quality cells.

Pei-Ying

I changed some things in the solver but didn't update the test case yet. For dynamic cases such as damBreak, change the Cpc coefficient to 0.5 in interFoam.C. Then the interface will move in a way that makes sense, but ironically there is too much smearing. At some point I'll also push a static testcase to the repository.

Cheers,

Anton

Hi, Anton,

Thanks for the reply!

I changed Cpc to 0.5 from 0.98, but, still got the same problem - interface did not move at all.

However, the code that I got several days ago seemed to get more reasonable results.

I am using OpenFOAM-2.1.x.

Pei-Ying

Hi All,

i ran the solver on of 2.0.0 - i had to edit the fvSolution and fvSchemes as Pei-Ying did.
When i run the dambreak scenario up to t = 5, and cpc=0.98 the water does not move, but the velocity profile looks kind of weird.

I'll start digging through your code..

EDIT: If I use CPC=0.5 it does not change anything at all...

Greetings

Attachment 12963

This is what it looks like for me (animated gif):
http://dl.dropbox.com/u/1261220/dambreak.gif

I'll push some more changes tomorrow.

- Anton

Hi Anton,

one question to your implementation: In the interfaceProperties.C on line 129, you define "factor". If I understood it right, the first part is the sqrt of Eq. 15, but what about the two other terms there? I couldn't find the connection to the paper...

regards

I have have made some changes in pimple calls in interFoam.C in order to compile it with OF2.1.0:

in line 93:
while (pimple.loop() || cycle<3)

in line 153:
for (int corr=0; corr<pimple.nCorrPISO(); corr++)

and also in pEqn.H, line 42:
mesh.solver(p.select(pimple.finalInnerIter()))

Also, I had to edit fvSolution and fvSchemes as Pei-Ying's post, and add the solvers for pcFinal and UFinal.

The damBreak testcase works fine.

Regards

Robert

In Eq. 15 you have sqrt(a*(1-a)), which is undefined for a=0 and a=1. I used the two extra terms to make sure the factor would become zero for 0,1, but it's probably not necessary.

I pushed updates to the code and the dambreak testcase, and added the static droplet testcase. The code now writes the surface tension force to file, and I believe it's sign in the static droplet testcase is wrong, but I'm not sure why. Also, the magnitude of spurious currents does not reduce to the levels of the reference paper in this test case...

Would any one of the OF-2.1.x users with a google account take the time to create a public clone for others to use?

- Anton

Ok thanks for the information..

maybe some critical works to the publication - it would be interesting what you think. To be very polemic, one could say, that he just interpolates and averages until the high velocities vanish?! I can't see a proof, that this is still "exact" what he does with the equations..

Has anyone tried (of those who get reasonable results) to set Cpc = 0 to get the original CSF formulation and compared it to the origianl interFoam solver?

Edit: With the new update, the dambreak works fine as well, so I'll give it a try

regards

I just realized sqrt(0) actually does work, so the factor here really does become unnecessary unless you want to reuse the result of 'w'. [edit] Problems also can arise when minimal unboundedness of alpha1 makes the term under the square root negative.

About your comments: Smoothing before calculating the curvature is known to be beneficial and has already been proposed by Brackbill. Separating the capillary pressure is only mentioned in the paper in conjunction with the filtering (which I haven't implemented yet), but just like the separation of rho*g I can imagine it to be beneficial for the numerics.

Hi,

I have created a clone of the code and modified the interFoam.C and pEqn.H to make it compatible with OF 2.1.0.

Also I have modified the dambreak testcase.

Robert