[rcastilla]

I have been working on capillary flows simulation with interFoam. When I make a simple simulation of the flow between two parallel vertical plates (2D), I got very strange results, due to an inadequate calculation of the contact angle correction (see the thread http://www.cfd-online.com/Forums/ope...tml#post264062 )

In the interfaceProperties library the contact angle correction is performed by modifying the direction of the unitary alpha gradient in the patch. It doesn't modify the alpha distribution itself, but only the vector that describes its gradient direction.

Then, in the interFoam solver, in the pEqn.H file, the contribution of the surface tension force is calculated as

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

The value of sigmaK is correct (it has been corrected by the interface object), but it is pointing in the wrong direction in the boundary cells, since it is calculated with fvc:snGrad(alpha1), which has not been corrected. The effect of that is huge when the flow is dominated by wall adhesion.

I suggest the following expression:

fvc::interpolate(interface.sigmaK()*mag(fvc::grad( alpha1)))*interface.nHatf()

It uses the gradient of alpha for the calculation of the force, but the direction is given by the surfaceScalarField nHatf, which has been corrected in the boundaries by the interface object.

I am testing it, and it seems to work. Simulations are very slow...

Perhaps there is a more elegant way of doing that.

Hope it will be useful.

Robert

[richpaj]

Surely though you would only want to apply the new formulation on the boundary wall and not everywhere else?

But otherwise, yes, it does look persuasive for a correction to the boundary values of the curvature force. Presumably, and probably more importantly, you would also want to correct the
boundary flux terms in the pEqn cf.

phi = phiU +
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(a lpha)
- ghf*fvc::snGrad(rho)
)*rUAf*mesh.magSf();

Nice detective work, well done!

[rcastilla]

In fact, the correction of this formulation is only made in the alphaContactAngle boundaries. I think that the numerical calculation of the flux is identical for the rest of the domain.

For the second comment, I am afraid that I have not explained well. It is the flux term that I have corrected. The original pEqn in my distribution is

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

and I have modified it as

phi = phiU +
(
fvc::interpolate(interface.sigmaK()*mag(fvc::grad( alpha1)))*interface.nHatf()
+ fvc::interpolate(rho)*(g & mesh.Sf())
)*rUAf;


Perhaps it should be also corrected in the UEqn file, but it only affects when the momentumPredictor switch is on, which is not my case so far.

Thank you, Richard, for your comments.

Robert

[henry]

fvc::interpolate(mag(fvc::grad( alpha1)))*interface.nHatf()

is not equivalent to

fvc::snGrad(alpha1)*mesh.magSf()

even in the bulk and the difference is important for numerical stability. Nevertheless you make a valid point concerning the direction implied by fvc::snGrad(alpha1) in the near-wall half-cell and a correction consistent with that applied to interface.nHatf() on the boundary should also be applied. Ideally this should be done through a specification of the boundary gradient of alpha1 but so far I have not found a completely reliable approach for this. I will reconsider the issue ....

H

[rcastilla]

Henry,

you are the expert. I am just now beginning to understand the FOAM language.

After your comment, may I suggest some kind of correction factor, calculated in interfaceProperties, that should be applied to fvc::snGrad in the dynamic equations?
This correction factor can be calculated with the old and new values of nHatfv. I am working in that, but I have not yet an adequate formulation.

With kind regards
Robert

[henry]

If you could work out what such a correction should be then it would be possible to work out what the gradient BC on alpha1 should be and then the snGrad(alpha1) would be correct wherever it is used. However, note that the cell gradient of alpha1 used to nHatfv uses the boundary value which is currently calculated using a zeroGradient condition which should also be corrected.

H

[rcastilla]

I don't understand completely your last comment. Also the snGrad is calculated with the boundary values of alpha1. What is the difference between the surface normal gradient calculated with snGrad and with grad(alpha1) and nHatfv?

You were right concerning to the stability of my previous formulation. It worked quite right for a while and then it began to show an strange behaviour. I will try with the new formulation correcting the value of snGrad.

Regards

Robert

[rcastilla]

A correcting factor for snGrad(alpha1) doesn't work, since, in the boundary, snGrad(alpha1) is identicaly null.

I am testing the following formulation:

in alphaEqnSubCycle.H, just after the correction of the contact angle:

volVectorField gradAlpha = fvc::grad(alpha1);

surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha);
surfaceScalarField magGradAlphaf = mag(gradAlphaf);



surfaceScalarField snGradAlphaS = fvc::snGrad(alpha1)*mesh.magSf();

surfaceScalarField nHatf=interface.nHatf();

//correct the value of snGradAlpha in the boundaries
forAll(snGradAlphaS.boundaryField(),patchi)
{
snGradAlphaS.boundaryField()[patchi]=magGradAlphaf.boundaryField()[patchi]*nHatf.boundaryField()[patchi];
}

and, then, in pEqn.H:

phi = phiU +
(
fvc::interpolate(interface.sigmaK())*snGradAlphaS
+ fvc::interpolate(rho)*(g & mesh.Sf())
)*rUAf;

That corrects the value of snGrad(alpha1) only in the boundaries. I have not yet results, but... makes it sense?

[rcastilla]

It doesn't work. The dynamics is not correct. The water column should rise, but it decreases. I have checked the value of snGradAlphaS in the boundaries and it is non-null. I cannot realize why it is not working properly... I keep thinking.

[henry]

I have just pushed a development to the curvature calculation and alphaContactAngle BCs into OpenFOAM-1.7.x which may help you. See the git log for details.

H

[rcastilla]

Henry,

I think I was too optimistic with the word "solved" in the thread's title. It looks like this problem is far from being solved.

I tried your solution, with the "limit none" keyword in the alphaContactAngle. The value of alpha1 in the wall went up to 1.2. But the results were not good. The water column still decrease when it should rise. I can try also with "limit alpha", but I don't think this is the problem.

I am focusing now with the compression procedure of the interface. It uses the value of nHatf, which is still not corrected in the cells next to the walls. I tried to make the curvature correction before the compression procedure, but it does not enhance the results. When there is not compression (cAlpha=0) the results are much better.


Robert

[henry]

To avoid unboundedness of alpha1 on the boundary I introduced optional limiting but this reduces the formal accuracy of the BC. If the unboundedness is causing a problem use "limit gradient" or "limit alpha".

In my capillary-rise tests the new BC improves the interface profile in the near-wall two cells and the accuracy of the rise height slightly (it was already fairly accurate). Note that you MUST use the correct corresponding BC for p_rgh on the walls for which the contact-angle BC is applied to alpha1 otherwise you will get a spurious flux through the wall.

H

[holger_marschall]

Originally Posted by henry

In my capillary-rise tests the new BC improves the interface profile in the near-wall two cells and the accuracy of the rise height slightly (it was already fairly accurate).

Hi Henry,

would you share these test cases on capillary-rise?

best regards,

[rcastilla]

Henry,

could you at least give some brief description of your test cases?

In my case there are 2D simulations of two vertical parallel plates, with a gap of 3x10^-4 m. Sigma is 0.072 N/m and theta=22 degrees. In order to keep the capillary column low and save computational space, the gravity is 10 times higher than usual, i.e. 98.1 m/s^2. With this configuration, analytical computations say that liquid should rise up to more than 4.5 mm. With an initial level of 4 mm and 64 cells between plates, I usually get a fall of the column, except when the interface compression is disable.

If you want I can send to you the cases files by e-mail.

With kind regards

Robert

[henry]

Take a look at the capillaryRise tutorial case in 1.7.1/x:

OpenFOAM-1.7.x/tutorials/multiphase/interFoam/laminar/capilaryRise

(please excuse the typo in the case name)

H

[rcastilla]

Hi, Henry,

I have run the simulation. Analytical calculations gives a result of average alpha = 0.5 (almost exactly). But in the tutorial the column rise up to less than 0.47. If interface compression is disable, it increases until 0.5, which is ok, but the interface is not sharp.

Have you any idea why the interface compression affects in such a way the wall adhesion simulation?

Robert

[duongquaphim]

Hi Henry,

I did use the new interfaceProperties library your posted in OF1.7.x in OF1.6 and successfully tested in a single processor. The contact angle closed to the wall looks much better.

However, I got an error of reading libinterfaceProperties.so when running the case in parallel:

[0] #0 Foam::error:rintStack(Foam::Ostream&) in "/fhome/hoang/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libOpenFOAM.so"
[0] #1 Foam::sigFpe::sigFpeHandler(int) in "/fhome/hoang/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libOpenFOAM.so"
[0] #2 ?? in "/lib64/libc.so.6"
[0] #3 Foam::tmp<Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> > Foam::mag<Foam::Vector<double>, Foam::fvsPatchField, Foam::surfaceMesh>(Foam::GeometricField<Foam::Vect or<double>, Foam::fvsPatchField, Foam::surfaceMesh> const&) in "/fhome/hoang/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libinterfaceProperties.so"
[0] #4 Foam::interfaceProperties::calculateK() in "/fhome/hoang/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libinterfaceProperties.so"
[0] #5 Foam::interfaceProperties::interfaceProperties(Foa m::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::IOdictionary const&) in "/fhome/hoang/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libinterfaceProperties.so"
[0] #6 main in "/fhome/hoang/OpenFOAM/OpenFOAM-1.6/applications/bin/linux64GccDPOpt/interFoam"
[0] #7 __libc_start_main in "/lib64/libc.so.6"
[0] #8 _start at /usr/src/packages/BUILD/glibc-2.9/csu/../sysdeps/x86_64/elf/start.S:116
[node16:20491] *** Process received signal ***
[node16:20491] Signal: Floating point exception (8)
[node16:20491] Signal code: (-6)
[node16:20491] Failing at address: 0x4320000500b
[node16:20491] [ 0] /lib64/libc.so.6 [0x2aadc03976e0]
[node16:20491] [ 1] /lib64/libc.so.6(gsignal+0x35) [0x2aadc0397645]
[node16:20491] [ 2] /lib64/libc.so.6 [0x2aadc03976e0]
[node16:20491] [ 3] /fhome/hoang/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libinterfaceProperties.so(_ZN4Foam3magINS_6VectorI dEENS_13fvsPatchFieldENS_11surfaceMeshEEENS_3tmpIN S_14GeometricFieldIdT0_T1_EEEERKNS6_IT_S7_S8_EE+0x 1d4) [0x2aadbd875c94]
[node16:20491] [ 4] /fhome/hoang/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libinterfaceProperties.so(_ZN4Foam19interfacePrope rties10calculateKEv+0x278) [0x2aadbd853d98]
[node16:20491] [ 5] /fhome/hoang/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libinterfaceProperties.so(_ZN4Foam19interfacePrope rtiesC1ERKNS_14GeometricFieldIdNS_12fvPatchFieldEN S_7volMeshEEERKNS1_INS_6VectorIdEES2_S3_EERKNS_12I OdictionaryE+0x635) [0x2aadbd854a05]
[node16:20491] [ 6] interFoam [0x4209a3]
[node16:20491] [ 7] /lib64/libc.so.6(__libc_start_main+0xe6) [0x2aadc0383586]
[node16:20491] [ 8] interFoam [0x41f6c9]
[node16:20491] *** End of error message ***
--------------------------------------------------------------------------
mpiexec noticed that process rank 0 with PID 20491 on node node16 exited on signal 8 (Floating point exception).
--------------------------------------------------------------------------

Have you ever tested this new implementation in parallel? Or can you please let me know how to fix this error?

Thank you.

Duong

[henry]

Thanks for the bug-report, this is now fixed in OpenFOAM-1.7.x

H

[alfa_8C]

Hy,
I'm using interFoam in a relatively simple setup. I have a tank filled with water. The top of the tank is open to atmosphere. Close to the bottom the tank has an outlet into a large resevoir of water. The free surface level of the reservoir is at -1m. As far as I know the BC for the pressure at the outlet has to be set to 0, as no hydrostatic pressure has to be considered. BUT!!! How can I make the code understand, that my free surface level outside the domain is at -1m? This must be given for sure but I don't know where.
Could anybody give me a hint?
Thanks in advance, Tony

[Andrea_85]

Hi all,

I've just updated my version of OF to 1.7.1. i'm still getting an error when the simulation is launched in parallel. (no error for single processor). What i have to do to fix the problem?

[1] #1 Foam::sigFpe::sigFpeHandler(int) in "/opt/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libOpenFOAM.so"
[2] #1 Foam::sigFpe::sigFpeHandler(int) in "/opt/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libOpenFOAM.so"
[1] #2 in "/opt/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libOpenFOAM.so"
[2] #2 ?? in "/lib64/libc.so.6"
[1] #3 void Foam::mag<Foam::Vector<double>, Foam::fvsPatchField, Foam::surfaceMesh>(Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh>&, Foam::GeometricField<Foam::Vector<double>, Foam::fvsPatchField, Foam::surfaceMesh> const&)?? in "/opt/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libinterfaceProperties.so"
[1] #4 Foam::tmp<Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> > Foam::mag<Foam::Vector<double>, Foam::fvsPatchField, Foam::surfaceMesh>(Foam::GeometricField<Foam::Vect or<double>, Foam::fvsPatchField, Foam::surfaceMesh> const&) in "/lib64/libc.so.6"
[2] #3 void Foam::mag<Foam::Vector<double>, Foam::fvsPatchField, Foam::surfaceMesh>(Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh>&, Foam::GeometricField<Foam::Vector<double>, Foam::fvsPatchField, Foam::surfaceMesh> const&) in "/opt/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libinterfaceProperties.so"
[1] #5 Foam::interfaceProperties::calculateK() in "/opt/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libinterfaceProperties.so"
[2] #4 Foam::tmp<Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> > Foam::mag<Foam::Vector<double>, Foam::fvsPatchField, Foam::surfaceMesh>(Foam::GeometricField<Foam::Vect or<double>, Foam::fvsPatchField, Foam::surfaceMesh> const&) in "/opt/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libinterfaceProperties.so"
[1] #6 Foam::interfaceProperties::interfaceProperties(Foa m::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::IOdictionary const&) in "/opt/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libinterfaceProperties.so"
[2] #5 Foam::interfaceProperties::calculateK() in "/opt/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libinterfaceProperties.so"
[1] #7 in "/opt/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libinterfaceProperties.so"
[2] #6 Foam::interfaceProperties::interfaceProperties(Foa m::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::IOdictionary const&) in "/opt/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libinterfaceProperties.so"
[2] #7 main in "/opt/OpenFOAM/OpenFOAM-1.7.1/applications/bin/linux64GccDPOpt/interFoam"
[1] #8 __libc_start_mainmain in "/lib64/libc.so.6"
[1] #9 in "/opt/OpenFOAM/OpenFOAM-1.7.1/applications/bin/linux64GccDPOpt/interFoam"
[2] #8 __libc_start_mainFoam::regIOobject::writeObject(Fo am::IOstream::streamFormat, Foam::IOstream::versionNumber, Foam::IOstream::compressionType) const in "/lib64/libc.so.6"
[2] #9 in "/opt/OpenFOAM/OpenFOAM-1.7.1/applications/bin/linux64GccDPOpt/interFoam"
[master:13645] *** Process received signal ***
[master:13645] Signal: Floating point exception (8)
[master:13645] Signal code: (-6)
[master:13645] Failing at address: 0x3f50000354d
[master:13645] [ 0] /lib64/libc.so.6 [0x2aee33d262d0]
[master:13645] [ 1] /lib64/libc.so.6(gsignal+0x35) [0x2aee33d26265]
[master:13645] [ 2] /lib64/libc.so.6 [0x2aee33d262d0]
[master:13645] [ 3] /opt/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libinterfaceProperties.so(_ZN4Foam3magINS_6VectorI dEENS_13fvsPatchFieldENS_11surfaceMeshEEEvRNS_14Ge ometricFieldIdT0_T1_EERKNS5_IT_S6_S7_EE+0x44) [0x2aee30e891e4]
[master:13645] [ 4] /opt/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libinterfaceProperties.so(_ZN4Foam3magINS_6VectorI dEENS_13fvsPatchFieldENS_11surfaceMeshEEENS_3tmpIN S_14GeometricFieldIdT0_T1_EEEERKNS6_IT_S7_S8_EE+0x 1a4) [0x2aee30e98134]
[master:13645] [ 5] /opt/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libinterfaceProperties.so(_ZN4Foam19interfacePrope rties10calculateKEv+0x2c5) [0x2aee30e74c75]
[master:13645] [ 6] /opt/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libinterfaceProperties.so(_ZN4Foam19interfacePrope rtiesC1ERKNS_14GeometricFieldIdNS_12fvPatchFieldEN S_7volMeshEEERKNS1_INS_6VectorIdEES2_S3_EERKNS_12I OdictionaryE+0x5e1) [0x2aee30e758e1]
[master:13645] [ 7] interFoam [0x420f2c]
[master:13645] [ 8] /lib64/libc.so.6(__libc_start_main+0xf4) [0x2aee33d13994]
[master:13645] [ 9] interFoam(_ZNK4Foam11regIOobject11writeObjectENS_8 IOstream12streamFormatENS1_13versionNumberENS1_15c ompressionTypeE+0xf1) [0x41dde9]
[master:13645] *** End of error message ***



Thanks

andrea