parasitic currents

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

 May 8, 2012, 08:58 #41 Senior Member   Andrea Ferrari Join Date: Dec 2010 Posts: 319 Rep Power: 15 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,377 Rep Power: 29 Code: hg clone https://code.google.com/r/robertocastillalopez-interfoamssf-210/ __________________ *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: 14
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:
 if (val[0] == "Time"): time = val[2] if (val[0] == "max(U):"): val[9] = float(val[9].replace("(", "")) print val[9] val[10] = float(val[10]) val[11] = float(val[11].replace(")", "")) #maxU = max(val[2], val[3], val[4]) maxU = math.sqrt(val[9]*val[9]+val[10]*val[10]+val[11]*val[11]) outfile.write(str(time) + "\t" + str(maxU) + "\n")
Greetings

 May 8, 2012, 09:32 #44 Senior Member     Anton Kidess Join Date: May 2009 Location: Germany Posts: 1,377 Rep Power: 29 The output of foamCalcEx looks like this: Code: Time = 0 Reading U Calculating maxU: (0 0 0) Of course you can let the python script parse the output of interFoamSSF itself, and then indeed you should search for max(U). - 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: 14 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: 319 Rep Power: 15 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,377
Rep Power: 29
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
Attached Images
 maxU.png (5.0 KB, 103 views)
__________________
*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: 319 Rep Power: 15 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 2-D. 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,377 Rep Power: 29 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: 319 Rep Power: 15 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:37 #52 Senior Member   Andrea Ferrari Join Date: Dec 2010 Posts: 319 Rep Power: 15 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,377
Rep Power: 29
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 OF-1.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:
 Originally Posted by Andrea_85 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=p-pc. Instead in your implementation there's p and there isn't explicitly any relationship between p, pd and pc. [...] best andrea
__________________
*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: 319 Rep Power: 15 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,377 Rep Power: 29 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: 319 Rep Power: 15 good job! you are right. andrea

 May 9, 2012, 14:37 #57 Senior Member   Andrea Ferrari Join Date: Dec 2010 Posts: 319 Rep Power: 15 Hi, i run again the static simulation e maximum velocities are now below 1e-5 after 0.002 s of simulations (with adjustable time step, it works!). it is not exactly the same as in the paper (they got 1e-8 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,377 Rep Power: 29 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 1e-5 at 5e-4s, but then diverge soon after (only in the liquid-gas 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: 319 Rep Power: 15 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 liquid-liquid case, i'll tell you tomorrow about the liquid-gas 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 1e-7 and then it rises up to 2e-6 (still lower than the paper, they have 40e-6s using 3 outerCorrector) while, using fixed time step, i have to use 1e-7 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,377
Rep Power: 29
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:
 Originally Posted by Andrea_85 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 liquid-liquid case, i'll tell you tomorrow about the liquid-gas 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 1e-7 and then it rises up to 2e-6 (still lower than the paper, they have 40e-6s using 3 outerCorrector) while, using fixed time step, i have to use 1e-7 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
__________________
*Spend as much time formulating your questions as you expect people to spend on their answer.

 Tags capillary flows, interfoam, parasitic currents