CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

"real" contact angle in the first cell layer

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   January 18, 2010, 10:56
Default "real" contact angle in the first cell layer
  #1
Member
 
Robert Castilla
Join Date: Apr 2009
Location: Spain
Posts: 80
Rep Power: 8
rcastilla is on a distinguished road
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.
Attached Images
File Type: jpg placas-3-0.1.jpg (18.3 KB, 123 views)
rcastilla is offline   Reply With Quote

Old   January 18, 2010, 14:49
Default
  #2
Senior Member
 
sega's Avatar
 
Sebastian Gatzka
Join Date: Mar 2009
Location: Frankfurt, Germany
Posts: 729
Rep Power: 11
sega is on a distinguished road
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.
__________________
Schrödingers wife: "What did you do to the cat? It's half dead!"
sega is offline   Reply With Quote

Old   January 18, 2010, 17:48
Default
  #3
Member
 
Robert Castilla
Join Date: Apr 2009
Location: Spain
Posts: 80
Rep Power: 8
rcastilla is on a distinguished road
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 is offline   Reply With Quote

Old   January 20, 2010, 09:50
Default Some results...
  #4
Member
 
Robert Castilla
Join Date: Apr 2009
Location: Spain
Posts: 80
Rep Power: 8
rcastilla is on a distinguished road
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?
Attached Images
File Type: jpg contact_angle.jpg (37.9 KB, 111 views)
rcastilla is offline   Reply With Quote

Old   January 21, 2010, 05:00
Default
  #5
Member
 
Richard Kenny
Join Date: Mar 2009
Posts: 59
Rep Power: 9
richpaj is on a distinguished road
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.
richpaj is offline   Reply With Quote

Old   January 21, 2010, 08:39
Default
  #6
Member
 
Robert Castilla
Join Date: Apr 2009
Location: Spain
Posts: 80
Rep Power: 8
rcastilla is on a distinguished road
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 parasitic currents), 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 is offline   Reply With Quote

Old   March 23, 2010, 16:08
Default
  #7
Member
 
Robert Castilla
Join Date: Apr 2009
Location: Spain
Posts: 80
Rep Power: 8
rcastilla is on a distinguished road
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.
rcastilla is offline   Reply With Quote

Old   March 24, 2010, 04:21
Default
  #8
ngj
Senior Member
 
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Deltares, Delft, The Netherlands
Posts: 1,619
Rep Power: 25
ngj will become famous soon enoughngj will become famous soon enough
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
ngj is offline   Reply With Quote

Old   March 25, 2010, 05:03
Default
  #9
Member
 
Robert Castilla
Join Date: Apr 2009
Location: Spain
Posts: 80
Rep Power: 8
rcastilla is on a distinguished road
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 is offline   Reply With Quote

Old   June 18, 2010, 06:06
Default
  #10
Member
 
Robert Castilla
Join Date: Apr 2009
Location: Spain
Posts: 80
Rep Power: 8
rcastilla is on a distinguished road
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 is offline   Reply With Quote

Old   June 22, 2010, 14:12
Default
  #11
Member
 
Robert Castilla
Join Date: Apr 2009
Location: Spain
Posts: 80
Rep Power: 8
rcastilla is on a distinguished road
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 is offline   Reply With Quote

Old   June 25, 2010, 06:55
Thumbs up
  #12
Member
 
Robert Castilla
Join Date: Apr 2009
Location: Spain
Posts: 80
Rep Power: 8
rcastilla is on a distinguished road
I think I have found the error in interFoam. I have open a new thread in OpenFoam/Bugs: solved: contact angle correction in interFoam

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
rcastilla is offline   Reply With Quote

Old   September 29, 2010, 12:39
Default Antidiffusion correction for sharp interface
  #13
Senior Member
 
Suresh kumar Kannan
Join Date: Mar 2009
Location: Luxembourg, Luxembourg, Luxembourg
Posts: 129
Rep Power: 8
kumar is on a distinguished road
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 is offline   Reply With Quote

Old   March 16, 2011, 05:04
Default Curvature models
  #14
Senior Member
 
Suresh kumar Kannan
Join Date: Mar 2009
Location: Luxembourg, Luxembourg, Luxembourg
Posts: 129
Rep Power: 8
kumar is on a distinguished road
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
kumar is offline   Reply With Quote

Old   March 17, 2011, 12:08
Default
  #15
ala
New Member
 
Alicja M
Join Date: Mar 2009
Location: Erlangen, DE
Posts: 26
Rep Power: 8
ala is on a distinguished road
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:
\nabla \alpha \approx \frac{1}{V} \int_V \nabla \alpha \approx \frac{1}{V} \sum_f {\alpha_f} \cdot S_f
where ()_f is interpolation to the surface (face of C.V.), V is your C.V.
_
Next: interpolation from cell center to the surface (\nabla \alpha)_f
_
Next: calculate free surface normal vector \hat{n}, the part at the face of C.V.
\hat{n}_f = \frac{\nabla \alpha_f}{|\nabla \alpha_f| + eps}
_
Next: calculate unit interface normal flux
nHatf = \phi_{\hat{n}} = \hat{n}_f \cdot S_f = \hat{n} \cdot n |S_f|
_
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

Last edited by ala; April 5, 2011 at 08:58.
ala is offline   Reply With Quote

Old   April 5, 2011, 09:09
Default
  #16
ala
New Member
 
Alicja M
Join Date: Mar 2009
Location: Erlangen, DE
Posts: 26
Rep Power: 8
ala is on a distinguished road
Hallo K.Suresh kumar,

in my opinion the curvature is estimated as mean value approximation, I mean with this following:
\kappa = \frac{1}{|V|} \int_V \nabla \cdot \hat{n} \: dV
          =  \frac{1}{|V|} \sum_f \int_{S_f} \hat{n} \cdot n \: dS
         \approx  \frac{1}{|V|} \sum_f |S_f| \hat{n} \cdot n

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

It is correct?

Greetings, Alicja
ala is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Dynamic contact angle rmousavibt Fluent UDF and Scheme Programming 10 March 7, 2014 08:00
Import netgen mesh to OpenFOAM hsieh Open Source Meshers: Gmsh, Netgen, CGNS, ... 32 September 13, 2011 05:50
Trimmed cell and embedded refinement mesh conversion issues michele OpenFOAM Other Meshers: ICEM, Star, Ansys, Pointwise, GridPro, Ansa, ... 2 July 15, 2005 04:15
Warning 097- AB CD-adapco 6 November 15, 2004 05:41
errors Fahad Main CFD Forum 0 March 23, 2004 14:20


All times are GMT -4. The time now is 18:00.