[rcastilla]

Hi,

I am trying to simulate capillary driven flows in micro-channels. When I saw some strange results, I made a simple test-bench: I simply simulated the capillary climbing of fluid between two parallel plates. Comparing with analytical results, the climbing of the simulated flow was lower, which implies an smaller contact angle than the one imposed in boundary conditions.

The contour of the interface is attached. It seems that up to the second cell layer, the angle is correct. But in the first layer the contact angle is almost double. I suspect that it is somehow related to spatial discretization, but I am not so skilled with C++ and Openfoam code to easily find and correct that.

I think that for most flows the difference is inappreciable. But for my 300microns channel, it is important.

I will thank any suggestion.

Regards

P.S. I have tried also with Fluent and the behaviour is more or less the same.

[sega]

How did you initialize the flow?
Does it correspond to the right contact angle?

Did you try different discretization?

Note that the contact angel which you specify in the boundary conditions is processed into a normal vector of the interface like this:

Normal vector of the interface = normal vector of the wall * cos(contact angle) + tangential vector of the wall * sin(contact angle).

So if you don't have "good" matching conditions at the second mesh layer (which may be governed by the flow itself and not the boundary condition) it may result is such a bad behaviour.

Tell us what happens if you spend more control volumes close to the wall.

[rcastilla]

Hi, Sebastian,

thank you for your answer.

Well, if I am not wrong, the definition of the normal vector is with the angle between interface and wall. That is:
normal vector to interface = normal vector to wall*sin(theta) + tangencial vector to wall*cos(theta)

Please, can you confirm if am wrong?

I can not understand how the conditions in the second layer can be bad, and how can it be related with the wrong behaviour in the first layer. The flow is governed by gravity and wall adhesion. In the interface shown, it was in equilibrium (i.e., there were no flow), and the weight of the fluid is such that the contact angle has to be (+ or -) double than the one specified in "0/alpha1". It seems that this value is used to calculate the correct shape of the interface starting in the second layer from the wall.

I can tell you that the behaviour changes with discretization, but I can not remember exactly how. I will check it tomorrow in the office and I will put here some results.

The initial field was with an horizontal interface, normal to walls.

Best regards

Robert

[rcastilla]

Here I post the results for three different mesh sizes. The labels are for the edgeGrading in blockMeshDict, with the same number of cells. The sizes of the cell close to the walls are:

D2/D1= 4 -> Delta_x=2.755e-3 mm
D2/D1 = 5 -> Delta_x=2.4e-3 mm
D2/D1 = 10 -> Delta_x=1.5e-3 mm

Te correct value of theta in the wall should be 22º. This angle is converging correctly towards the wall, except in the last node, when is jump to a higher value, which implies a lower wall adhesion force and, hence, a lower surface level, as shown in the left figure.

Has somebody any idea of the reason of this strange behaviour?

[richpaj]

This is my understanding of what is happening, hopefully not too confusingly explained.

The contact angle information (as specified in the gamma or alpha field) is incorporated into the estimation of curvature which is subsequently used to solve the U and p equations. The profile of gamma (or alpha if you're referring to OF16x) is thus derived from the behaviour of the dynamics i.e. via U when the gammaEqn is solved.

So, the above is fine up to the cells adjoining the wall.

At the wall we require the contact line to move (!) (hence the use of the
patchInterField() values for U which are cell values adjacent to the wall)
so that, for example, we can estimate a velocity-dependent contact-angle.

In this scheme, therefore, it appears that the value of gamma exactly on the wall is not
in fact used.
This possibly explains why "gammaContactAngleFvPatchScalarField" is derived
from a "zeroGradientFvPatchField" which forces the boundary value of gamma to be the same as that of the adjacent boundary cells. So the boundary value is somewhat arbitrarily (but reasonably) specified.

If your contact angles are correct up to the cells adjacent to the wall then that should be fine.
The relevant contact angle information is being incorporated into the dynamics despite the fact
it might not be aesthetically pleasing exactly on the wall itself.

Another aspect to "getting the correct contact angle" is the resolution required to ensure that we have decent values of U near to the wall.
This can be ascertained by examining the equivalent of a yPlus which should be ~1.

Good luck,

Richard K.

[rcastilla]

Thank you, Richard, for your explanation. It makes sense.

I am thinking that, in fact, there is no U, since the flow is in equilibrium, and the interface position is constant, but there are the parasitic currents (see the thread http://www.cfd-online.com/Forums/ope...-currents.html), and these could be the reason of this behaviour.

But the problem is that the dynamics (in fact, statics) is not correct. This is not only an aesthetic question. The position of the interface does not agree with the analytically calculated.

So, I should make a very thin layer adjacent to the wall I'll try, reducing the parasitic velocities. I will post the results.

[rcastilla]

This is the state of the problem:

• Two vertical plates separated by a gap of 300 microns
• Gravity multiplied by 10 (i.e. g=98.1 m/s^2) in order to reduce the computational domain.
• pressure zero in the bottom and the top, so that flow is driven only by gravity and wall adhesion
• contact angle = 22 degrees and surface tension coefficient = 0.072 N/m
• Analytical calculation gives a mean height of the water column of, approx, 4.5 mm in static equilibrium.
• Numerical results give smaller water columns, which are in agreement with a bigger contact angle, as shown in previous post.

I have tried with better space discretitation near the wall, but the results are even worse.

I have played then with different numerical parameters, and, if the interface compression factor (cAlpha) is 0, i.e. no compression at all, the contact angle seems to be correct, and the column height is about 4.3 mm, with is in good agreement with analytical results. The interface, though, looks very smeared.

That makes sense. I think that, since the interface compression is made with an articifial diffusion term with a velocity, which is obtained with phi:

surfaceScalarField phic = mag(phi/mesh.magSf());
phic = min(interface.cAlpha()*phic, max(phic));
surfaceScalarField phir = phic*interface.nHatf();

then, near the wall this "artifical" velocity vanishes, and this compression process distorts the interface and, hence, the contact angle.

What I don't understand is why this contact angle is not corrected by the instruction interface.correct() just after the solution of the alpha equation in alphaEqnSubCycle.H. This routine is made precisely to correct the contact angle in the first layer in the wall.

Any suggestion will be welcome.

[ngj]

Hi Robert

I does indeed sound like the VOF-method is effecting the results. I do have an alternative (much less than beautifully implemented) VOF-solver, which it could be fun to test. I do not know if it would work with contact angle, etc. However send me an email abd we might find a way to test your problem.

Best regards,

Niels

[rcastilla]

Actually, I don't think that the problem is in the VOF method, but only in the interface compression procedure.

I have found this reference:

K. K. So, X. Y. Hu and N. A. Adams, Anti-diffusion correction for sharp interface of two-phase incompressible flow, Open Source CFD International Conference 2009

that suggest an alternative method for the interface compression. Has anybody tested it with wall adhesion?

[rcastilla]

Hi,

after many tests, I concluded that the interface compression modifies, somehow, the interface angle in the first cell attached to the wall.

This modification depends on the mesh size. I have tested with different mesh sizes in the normal direction to the wall. The behavior can be very strange. Keeping constant the mesh except in the first cell, it seems to behave better than with a grading mesh (with the same size in the first cell!).

I have checked the interface compression code and the interface contact angle code, but I have no idea how the first can affect the second.

Can someone shed some light on this problem?

BTW, I have tested also the alternative compression algorithm by So et al. but the behaviour is more or less the same.

Robert

[rcastilla]

Finally, I think that the interface compression procedure is not the cause of the modification of the contact angle, but, as it modifies the gradient of alpha, the effect is different. In fact, it is worse when the gradient is larger, i.e., when there is an interface compression.

It seems that Kennys's post (#5) is the key. Indeed, it looks like the angle in the patch is not actually used, because of the zeroGradient boundary condition in alpha. The dynamics of the flow is controlled by the interface angle in the next face close to the wall. I think that a possible workaround could be to impose the contact angle in the cell next to the patch. I could do that with faceCells() if the contact angle were a constant value, but it is a field. Maybe with a loop over all the faces in the patch and then with owner(). But then I have the problem that in interfaceProperties the nHat is defined, and corrected, in the faces:

const fvMesh& mesh = alpha1_.mesh();
const surfaceVectorField& Sf = mesh.Sf();

// Cell gradient of alpha
volVectorField gradAlpha = fvc::grad(alpha1_);

// Interpolated face-gradient of alpha
surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha);
//gradAlphaf -=
// (mesh.Sf()/mesh.magSf())
// *(fvc::snGrad(alpha1_) - (mesh.Sf() & gradAlphaf)/mesh.magSf());

// Face unit interface normal
surfaceVectorField nHatfv = gradAlphaf/(mag(gradAlphaf) + deltaN_);
correctContactAngle(nHatfv.boundaryField());


It should be possible to change the value of nHatfv in the faces next to the patch, but I am not able to find a way to perform a loop over the faces of a cell.

I keep thinking, but any suggestion will be wellcome.

Robert

[rcastilla]

I think I have found the error in interFoam. I have open a new thread in OpenFoam/Bugs: http://www.cfd-online.com/Forums/ope...interfoam.html

Again the key was the post by Kenny. In despite of the shape of the interface in the wall, the force should be correctly incorporated into the dynamics. I am testing the corrected version and it seems to work , but it is still too soon. Simulations are very slow....

I hope it will be useful. Thank you for the help.

Robert

[kumar]

Hello rcastilla,
recently I came across the paper "Antidiffusion correction for sharp interface of two-phase incompressible flow"

I am using OF-1.6.x and I would like to test this antidiffusion method for my problem.
I read that you had already tried this antidffusion method.

Could you let me know the steps you did to include the antidiffusion equation into the interFoam solver.

Thankyou very much
regards
K.Suresh kumar

[kumar]

Hello,
I just got hold of the paper "A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework" Journal of Computational Physics 213 (2006) 141-173.

One of the major conclusions of the paper is that the origin of spurious currents, regardless of the surface tension model employed, is errors in curvature estimates.

They also explain two curvature estimation models:
Convolution technique and Height function technique.

Can anyone tell me in simple formulation how the curvature is estimated in interFoam ?

Thanks
regards
K.Suresh kumar

[ala]

Hallo,

take a look at this thread:
vof-method

I think, there is no convolution/smoothing of VOF-function alpha. Because of numeric diffusion, you have something like smooth VOF-function alpha.

In implementation interFoam, if you use compressionScheme, you are trying to keep your function alpha so smooth as possible.

The implementation of curvature term is in file interfaceProperties.C, see: calculateK().

I'm also interested to know estimation of curvature in interFoam.

______

So far, i can read from this file:
_
Firstly you calculates:

where is interpolation to the surface (face of C.V.), is your C.V.
_
Next: interpolation from cell center to the surface
_
Next: calculate free surface normal vector , the part at the face of C.V.

_
Next: calculate unit interface normal flux
nHatf =
_
Next: calculate divergence
K = -fvc::div(nHatf)
---> See doxygen .
_

Please, take a look at this, and if i'm wrong tell it to me.

Greetings, Alicja

[ala]

Hallo K.Suresh kumar,

in my opinion the curvature is estimated as mean value approximation, I mean with this following:


How is the capillary force term in UEqn.H and pEqn.H discretised, is an another story.

It is correct?

Greetings, Alicja