
[Sponsors] 
May 8, 2012, 08:58 

#41 
Senior Member
Andrea Ferrari
Join Date: Dec 2010
Posts: 314
Rep Power: 10 
Hi Robert ,
how the version that compiles with OF 2.1 can be downloaded? I hope to have time to test it as soon as possible! thanks andrea 

May 8, 2012, 09:01 

#42 
Senior Member
Anton Kidess
Join Date: May 2009
Location: Germany
Posts: 1,267
Rep Power: 23 
Code:
hg clone https://code.google.com/r/robertocastillalopezinterfoamssf210/
__________________
*On twitter @akidTwit *Spend as much time formulating your questions as you expect people to spend on their answer. 

May 8, 2012, 09:10 

#43  
Senior Member
Join Date: Nov 2010
Posts: 113
Rep Power: 9 
Hi Anton,
thanks for your opinion. Just one comment to your update: The maxU.py script searches for "Calculating" to extract the velocities.. but this only occurs for the gravity, or did i get it wrong? I let it search for max(U).. this is my edit: Quote:


May 8, 2012, 09:32 

#44 
Senior Member
Anton Kidess
Join Date: May 2009
Location: Germany
Posts: 1,267
Rep Power: 23 
The output of foamCalcEx looks like this:
Code:
Time = 0 Reading U Calculating maxU: (0 0 0)  Anton
__________________
*On twitter @akidTwit *Spend as much time formulating your questions as you expect people to spend on their answer. 

May 8, 2012, 09:33 

#45 
Senior Member
Join Date: Nov 2010
Posts: 113
Rep Power: 9 
ah ok, sorry, now i got how your script was supposed to run
Thanks! 

May 8, 2012, 11:51 

#46 
Senior Member
Andrea Ferrari
Join Date: Dec 2010
Posts: 314
Rep Power: 10 
Hi,
just one question: why the surface force are first calculated at the face centres and then averaged back to the cell centres before solving pc equation in your implementation? I think the idea behind the balance force algorithm is to calculate both capillary pressure and surface force at same place (cell centres or face centres) to reduce errors due to interpolation. So, is there a reason why you write: fvScalarMatrix pcEqn ( fvm::laplacian(pc) == fvc::div(fc) ); instead of vScalarMatrix pcEqn ( fvm::laplacian(pc) == fvc::div(fcf*mesh.magSf()) ); (and same for uEqn and pEqn). This is for example what is actually done in the standard interFoam solver. The surface forces are directly calculated at the face centres and not interpolated from cell centres values. andrea 

May 8, 2012, 12:53 

#47 
Senior Member
Anton Kidess
Join Date: May 2009
Location: Germany
Posts: 1,267
Rep Power: 23 
Excellent question Andrea! In fact, I've been testing that this afternoon and the capillary pressure field looks a lot better. I'll push a fix soon. Unfortunately the spurious currents are still larger than in the reference (O(0.01)), so keep your eyes open for more bugs
__________________
*On twitter @akidTwit *Spend as much time formulating your questions as you expect people to spend on their answer. 

May 8, 2012, 15:44 

#48 
Senior Member
Andrea Ferrari
Join Date: Dec 2010
Posts: 314
Rep Power: 10 
Hi,
i spent some time to test the drop relaxation and the dambreak case. I got very strange results in the static case where the interface starts to deform in an anomalous way and it is finally completely destroyed. Do you get the same results? or am i doing something wrong? because in this case i think it is more important understand the reason of that strange behavior than worrying about spurious currents...the simulation looks very bad. I also think that the dynamic case works fine only because the surface tension is much less important than in the static case, tomorrow i'll tryit again with dimension of mm or \mu m and i'm quite sure to get more or less the same strange behavior of the interface. I changed a little bit the test case, i made it 2D. It looks like there's a wrong sign somewhere, because of on the corners of the square the surface force should points inward (correct me if i am wrong) but it points outward and the smearing of the droplet starts actually from those corners. Do you agree on that? best andrea 

May 8, 2012, 16:14 

#49 
Senior Member
Anton Kidess
Join Date: May 2009
Location: Germany
Posts: 1,267
Rep Power: 23 
I agree there seems to be a sign error (see post #37). The image I posted earlier was for results where I changed the sign on deltasf to fix the issue. Then you get the proper shape for the droplet and the surface tension force vectors make sense. Why I have to add a minus sign I don't understand yet, and thus I haven't pushed that change to the repository.
 Anton
__________________
*On twitter @akidTwit *Spend as much time formulating your questions as you expect people to spend on their answer. 

May 8, 2012, 16:26 

#50 
Senior Member
Andrea Ferrari
Join Date: Dec 2010
Posts: 314
Rep Power: 10 
ok,
i missed the sign change you mentioned in your post. the equations look fine, may be the wrong minus comes from somewhere else. i have to write down the equations and think about that. best andrea 

May 9, 2012, 08:24 

#51 
Senior Member
Andrea Ferrari
Join Date: Dec 2010
Posts: 314
Rep Power: 10 
Hi again,
i have two comments 1) in createFields.H maybe it is better to replace surfaceScalarField fcf_old = interface.sigma()*fvc::snGrad(alpha1)*fvc::interpo late(interface.K()); with surfaceScalarField fcf_old = interface.sigma()*fvc::snGrad(alpha1)*interface.Kf (); to be consistent with the rest. 2)I am not sure about the implementation of the equations in terms of p or pd. In the paper there's pd, which is pd=ppc. Instead in your implementation there's p and there isn't explicitly any relationship between p, pd and pc. In particulary when you define the face fluxes you have: phi_star=(H/A)_f + (rho g/A)_f + (fc/A)_f  (1/A*grad pc)_f and phi = phi_star (1/A*grad p)_f then you solve for p by applying the divergence at the last equation and using the continuity. So i think that the final velocity: u=H/A + (rho*g)/A + (fc/A) (1/A*grad pc) (1/A)*grad p has one more term (grad pc) that there should not be there (grad pc is an artificial additional term, its contribution should already be in p, which is the total pressure). Solving for pd makes more sense in my opinion. You define pd in createFields as pd=ppc then you will have phi_star=(H/A)_f + (rho g/A)_f + (fc/A)_f  (grad pc/A)_f then you solve for pd: \nabla[1/A(grad pd)]=\nabla(phi_star) and update the face fluxes and total pressure p=pc+pd; phi=phi_star (1/A)*grad pd=phi_star 1/A*grad p + 1/A*grad pc and the contribution of grad pc should disappear from the velocity. Don't know if it makes sense or not, i spent all this morning thinking about that. what are your comments? best andrea 

May 9, 2012, 08:37 

#52 
Senior Member
Andrea Ferrari
Join Date: Dec 2010
Posts: 314
Rep Power: 10 
ok,
i probably just realized that your "p" is what is called "pd" in the paper. you simply do not define any total pressure. right? andrea 

May 9, 2012, 09:31 

#53  
Senior Member
Anton Kidess
Join Date: May 2009
Location: Germany
Posts: 1,267
Rep Power: 23 
1) I agree that it's better to be consistent, so I'll update that. The effect on the solution should be negligible though.
2) I'm guessing the authors of the paper started working with OF1.5 and thus have the pd term. In case of the interFoamSSF solver indeed I have been using 'p' like 'pd' and do not determine the total pressure. I'm sorry that this caused some confusion. Do you agree that the terms are the way they should be taking the different naming into account?  Anton Quote:
__________________
*On twitter @akidTwit *Spend as much time formulating your questions as you expect people to spend on their answer. Last edited by akidess; May 9, 2012 at 09:33. Reason: shortened quote 

May 9, 2012, 10:04 

#54 
Senior Member
Andrea Ferrari
Join Date: Dec 2010
Posts: 314
Rep Power: 10 
Hi Anton
if "pd" is what you called "p", then i think it is correct. One can simply define the total pressure like ptot=p+pc. I still thinking about the minus... best andrea 

May 9, 2012, 10:17 

#55 
Senior Member
Anton Kidess
Join Date: May 2009
Location: Germany
Posts: 1,267
Rep Power: 23 
I think I found it! Standard interFoam calculates the curvature as K_=fvc::div(nHatf_), whereas the paper doesn't have the minus sign (eq. 14)!
__________________
*On twitter @akidTwit *Spend as much time formulating your questions as you expect people to spend on their answer. 

May 9, 2012, 10:33 

#56 
Senior Member
Andrea Ferrari
Join Date: Dec 2010
Posts: 314
Rep Power: 10 
good job! you are right.
andrea 

May 9, 2012, 14:37 

#57 
Senior Member
Andrea Ferrari
Join Date: Dec 2010
Posts: 314
Rep Power: 10 
Hi,
i run again the static simulation e maximum velocities are now below 1e5 after 0.002 s of simulations (with adjustable time step, it works!). it is not exactly the same as in the paper (they got 1e8 after 0.001 s of sim) but it seems we are in the right direction. the difference might be due to filtering? I also noticed that they used the kernel to smooth the curvature also to smooth the gradient of the indicator function(what is called gradAlpha is your implementation), which is then used to get the interface normal (pag 7 of the paper). don't know if it might be important. do you get same results? best andrea 

May 9, 2012, 15:19 

#58 
Senior Member
Anton Kidess
Join Date: May 2009
Location: Germany
Posts: 1,267
Rep Power: 23 
If filtering is applied we should get the curve labeled FSF in the paper I believe, so that can't be it.
I also noticed the sentence on the smoothing of the VOF gradient, but I'm not sure what to do with it. It's a bit ambiguous in my opinion on how the kernel is applied  for every iteration of eq 12, or just once at the end? I'll test that tomorrow. About the residuals  I have them drop to 1e5 at 5e4s, but then diverge soon after (only in the liquidgas case). Do you see the same behavior? The problem I had with adjustable runTime was that OpenFoam wouldn't exactly half the time step for the subcycle, but adapt it to fit the write schedule  which of course is bad if you're relying on it to return half the flux. Are you sure it works the way you think it does?
__________________
*On twitter @akidTwit *Spend as much time formulating your questions as you expect people to spend on their answer. 

May 9, 2012, 16:31 

#59 
Senior Member
Andrea Ferrari
Join Date: Dec 2010
Posts: 314
Rep Power: 10 
I agree that is a bit ambiguous, i would apply it once at the end so
1) smooth alpha 2 times with eq 12 2)calculate the gradient using alpha2s_ 3) apply eq 15 with "gradAlpha" insted of "k" and gradAlpha_s* instead of "ks*" At the moment i tested only the liquidliquid case, i'll tell you tomorrow about the liquidgas case. (with "residual" do you mean maximum velocity or simualtion residual?) about the time step...I do not know exactly what OF does to fit the adjustable time step to the "writing time step", does it perform a sort of interpolation between two times? I this case i think it is not correct to take half the flux, but should i have problem with conservation of mass or not? the mass is conserved and the sim seems stable, so how can i check if it works the way it is supposed to work? with adjustable time step it starts with 1e7 and then it rises up to 2e6 (still lower than the paper, they have 40e6s using 3 outerCorrector) while, using fixed time step, i have to use 1e7 s to avoid floating point errors in the square root calculation due to little unboundedness of alpha (may this can be fix in another way, using for example a limited alpha). best andrea 

May 10, 2012, 03:09 

#60  
Senior Member
Anton Kidess
Join Date: May 2009
Location: Germany
Posts: 1,267
Rep Power: 23 
With "residual" I meant the remaining mean maximum velocity.
About the time stepping  print the time step and current time after each call of setDeltaT and setTime (using runTime.time().value() and runTime.deltaTValue()). You should not see the values you would expect. What did you set for "writeControl"? Limiting alpha1 in the calculation of "factor" is an effective way to get rid of the negative square root.  Anton Quote:
__________________
*On twitter @akidTwit *Spend as much time formulating your questions as you expect people to spend on their answer. 

Tags 
capillary flows, interfoam, parasitic currents 
Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
how to monitor free surface elevation vs time in OF?  ozgur  OpenFOAM PostProcessing  56  September 14, 2015 08:11 
parasitic currents  PeiYing Hsieh  Main CFD Forum  0  January 13, 2009 20:58 
Parasitic currents reduction  hsieh  OpenFOAM Running, Solving & CFD  0  January 13, 2009 16:44 
Parasitic currents reduction  hsieh  OpenFOAM Running, Solving & CFD  0  January 13, 2009 16:37 
Modelling ocean currents of the past Earth  pgm  Main CFD Forum  3  March 2, 2005 09:45 