|
[Sponsors] |
solved: contact angle correction in interFoam |
|
LinkBack | Thread Tools | Search this Thread | Display Modes |
June 25, 2010, 04:14 |
solved: contact angle correction in interFoam
|
#1 |
Senior Member
Robert Castilla
Join Date: Apr 2009
Location: Spain
Posts: 109
Rep Power: 17 |
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 |
|
June 25, 2010, 08:50 |
|
#2 |
Member
Richard Kenny
Join Date: Mar 2009
Posts: 64
Rep Power: 18 |
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! |
|
June 25, 2010, 09:48 |
|
#3 |
Senior Member
Robert Castilla
Join Date: Apr 2009
Location: Spain
Posts: 109
Rep Power: 17 |
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 |
|
June 27, 2010, 11:32 |
|
#4 |
Senior Member
Join Date: Mar 2009
Posts: 854
Rep Power: 22 |
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 |
|
June 28, 2010, 10:07 |
|
#5 |
Senior Member
Robert Castilla
Join Date: Apr 2009
Location: Spain
Posts: 109
Rep Power: 17 |
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 |
|
June 28, 2010, 10:20 |
|
#6 |
Senior Member
Join Date: Mar 2009
Posts: 854
Rep Power: 22 |
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 |
|
June 28, 2010, 11:12 |
|
#7 |
Senior Member
Robert Castilla
Join Date: Apr 2009
Location: Spain
Posts: 109
Rep Power: 17 |
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 |
|
June 29, 2010, 07:52 |
|
#8 |
Senior Member
Robert Castilla
Join Date: Apr 2009
Location: Spain
Posts: 109
Rep Power: 17 |
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? |
|
June 29, 2010, 18:07 |
|
#9 |
Senior Member
Robert Castilla
Join Date: Apr 2009
Location: Spain
Posts: 109
Rep Power: 17 |
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.
|
|
August 6, 2010, 07:55 |
|
#10 |
Senior Member
Join Date: Mar 2009
Posts: 854
Rep Power: 22 |
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 |
|
August 11, 2010, 06:35 |
|
#11 |
Senior Member
Robert Castilla
Join Date: Apr 2009
Location: Spain
Posts: 109
Rep Power: 17 |
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 |
|
August 14, 2010, 13:52 |
|
#12 |
Senior Member
Join Date: Mar 2009
Posts: 854
Rep Power: 22 |
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 |
|
August 15, 2010, 04:27 |
|
#13 | |
Senior Member
|
Quote:
would you share these test cases on capillary-rise? best regards,
__________________
Holger Marschall web: http://www.holger-marschall.info mail: holgermarschall@yahoo.de |
||
September 1, 2010, 06:06 |
|
#14 |
Senior Member
Robert Castilla
Join Date: Apr 2009
Location: Spain
Posts: 109
Rep Power: 17 |
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 |
|
September 1, 2010, 06:20 |
|
#15 |
Senior Member
Join Date: Mar 2009
Posts: 854
Rep Power: 22 |
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 |
|
September 2, 2010, 08:45 |
|
#16 |
Senior Member
Robert Castilla
Join Date: Apr 2009
Location: Spain
Posts: 109
Rep Power: 17 |
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 |
|
September 7, 2010, 12:51 |
|
#17 |
Member
|
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 |
|
September 7, 2010, 18:13 |
|
#18 |
Senior Member
Join Date: Mar 2009
Posts: 854
Rep Power: 22 |
Thanks for the bug-report, this is now fixed in OpenFOAM-1.7.x
H |
|
April 19, 2011, 05:22 |
|
#19 |
Member
Antonio Liggieri
Join Date: Aug 2010
Posts: 76
Rep Power: 15 |
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 |
|
April 19, 2011, 11:49 |
|
#20 |
Senior Member
Andrea Ferrari
Join Date: Dec 2010
Posts: 319
Rep Power: 16 |
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 |
|
Thread Tools | Search this Thread |
Display Modes | |
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
Dynamic contact angle | rmousavibt | Fluent UDF and Scheme Programming | 12 | October 31, 2021 22:38 |
Slug Flow, interFoam, problems with Contact Angle | PrzemekPL | OpenFOAM Running, Solving & CFD | 13 | February 18, 2014 22:10 |
[Netgen] Import netgen mesh to OpenFOAM | hsieh | OpenFOAM Meshing & Mesh Conversion | 32 | September 13, 2011 05:50 |
Implemening Slip bc and Dynamic contact angle in interFoam | asaha | OpenFOAM Running, Solving & CFD | 2 | July 29, 2010 15:49 |
Theoretical background of formula for dynamic contact angle in interfoam | sebastian_vogl | OpenFOAM Running, Solving & CFD | 3 | June 22, 2009 12:25 |