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

parasitic currents

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

Like Tree22Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   May 8, 2012, 08:58
Default
  #41
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 275
Rep Power: 7
Andrea_85 is on a distinguished road
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
Andrea_85 is offline   Reply With Quote

Old   May 8, 2012, 09:01
Default
  #42
Senior Member
 
akidess's Avatar
 
Anton Kidess
Join Date: May 2009
Location: Delft, Netherlands
Posts: 919
Rep Power: 17
akidess will become famous soon enough
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.
*Help define the OpenFOAM stackexchange Q&A site: http://area51.stackexchange.com/prop...oam-technology
akidess is offline   Reply With Quote

Old   May 8, 2012, 09:10
Default
  #43
Senior Member
 
Join Date: Nov 2010
Posts: 113
Rep Power: 6
lindstroem is on a distinguished road
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
lindstroem is offline   Reply With Quote

Old   May 8, 2012, 09:32
Default
  #44
Senior Member
 
akidess's Avatar
 
Anton Kidess
Join Date: May 2009
Location: Delft, Netherlands
Posts: 919
Rep Power: 17
akidess will become famous soon enough
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.
*Help define the OpenFOAM stackexchange Q&A site: http://area51.stackexchange.com/prop...oam-technology
akidess is offline   Reply With Quote

Old   May 8, 2012, 09:33
Default
  #45
Senior Member
 
Join Date: Nov 2010
Posts: 113
Rep Power: 6
lindstroem is on a distinguished road
ah ok, sorry, now i got how your script was supposed to run

Thanks!
lindstroem is offline   Reply With Quote

Old   May 8, 2012, 11:51
Default
  #46
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 275
Rep Power: 7
Andrea_85 is on a distinguished road
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
Andrea_85 is offline   Reply With Quote

Old   May 8, 2012, 12:53
Default
  #47
Senior Member
 
akidess's Avatar
 
Anton Kidess
Join Date: May 2009
Location: Delft, Netherlands
Posts: 919
Rep Power: 17
akidess will become famous soon enough
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
File Type: png maxU.png (5.0 KB, 64 views)
__________________
*On twitter @akidTwit
*Spend as much time formulating your questions as you expect people to spend on their answer.
*Help define the OpenFOAM stackexchange Q&A site: http://area51.stackexchange.com/prop...oam-technology
akidess is offline   Reply With Quote

Old   May 8, 2012, 15:44
Default
  #48
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 275
Rep Power: 7
Andrea_85 is on a distinguished road
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
Andrea_85 is offline   Reply With Quote

Old   May 8, 2012, 16:14
Default
  #49
Senior Member
 
akidess's Avatar
 
Anton Kidess
Join Date: May 2009
Location: Delft, Netherlands
Posts: 919
Rep Power: 17
akidess will become famous soon enough
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.
*Help define the OpenFOAM stackexchange Q&A site: http://area51.stackexchange.com/prop...oam-technology
akidess is offline   Reply With Quote

Old   May 8, 2012, 16:26
Default
  #50
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 275
Rep Power: 7
Andrea_85 is on a distinguished road
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
Andrea_85 is offline   Reply With Quote

Old   May 9, 2012, 08:24
Default
  #51
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 275
Rep Power: 7
Andrea_85 is on a distinguished road
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. 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=p-pc

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

Old   May 9, 2012, 08:37
Default
  #52
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 275
Rep Power: 7
Andrea_85 is on a distinguished road
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
Andrea_85 is offline   Reply With Quote

Old   May 9, 2012, 09:31
Default
  #53
Senior Member
 
akidess's Avatar
 
Anton Kidess
Join Date: May 2009
Location: Delft, Netherlands
Posts: 919
Rep Power: 17
akidess will become famous soon enough
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 View Post
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
__________________
*On twitter @akidTwit
*Spend as much time formulating your questions as you expect people to spend on their answer.
*Help define the OpenFOAM stackexchange Q&A site: http://area51.stackexchange.com/prop...oam-technology

Last edited by akidess; May 9, 2012 at 09:33. Reason: shortened quote
akidess is offline   Reply With Quote

Old   May 9, 2012, 10:04
Default
  #54
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 275
Rep Power: 7
Andrea_85 is on a distinguished road
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
Andrea_85 is offline   Reply With Quote

Old   May 9, 2012, 10:17
Default
  #55
Senior Member
 
akidess's Avatar
 
Anton Kidess
Join Date: May 2009
Location: Delft, Netherlands
Posts: 919
Rep Power: 17
akidess will become famous soon enough
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.
*Help define the OpenFOAM stackexchange Q&A site: http://area51.stackexchange.com/prop...oam-technology
akidess is offline   Reply With Quote

Old   May 9, 2012, 10:33
Default
  #56
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 275
Rep Power: 7
Andrea_85 is on a distinguished road
good job! you are right.

andrea
Andrea_85 is offline   Reply With Quote

Old   May 9, 2012, 14:37
Default
  #57
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 275
Rep Power: 7
Andrea_85 is on a distinguished road
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
Andrea_85 is offline   Reply With Quote

Old   May 9, 2012, 15:19
Default
  #58
Senior Member
 
akidess's Avatar
 
Anton Kidess
Join Date: May 2009
Location: Delft, Netherlands
Posts: 919
Rep Power: 17
akidess will become famous soon enough
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.
*Help define the OpenFOAM stackexchange Q&A site: http://area51.stackexchange.com/prop...oam-technology
akidess is offline   Reply With Quote

Old   May 9, 2012, 16:31
Default
  #59
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 275
Rep Power: 7
Andrea_85 is on a distinguished road
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
Andrea_85 is offline   Reply With Quote

Old   May 10, 2012, 03:09
Default
  #60
Senior Member
 
akidess's Avatar
 
Anton Kidess
Join Date: May 2009
Location: Delft, Netherlands
Posts: 919
Rep Power: 17
akidess will become famous soon enough
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 View Post
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
__________________
*On twitter @akidTwit
*Spend as much time formulating your questions as you expect people to spend on their answer.
*Help define the OpenFOAM stackexchange Q&A site: http://area51.stackexchange.com/prop...oam-technology
akidess is offline   Reply With Quote

Reply

Tags
capillary flows, interfoam, parasitic currents

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
how to monitor free surface elevation vs time in OF? ozgur OpenFOAM Post-Processing 55 October 31, 2013 10:33
parasitic currents Pei-Ying 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


All times are GMT -4. The time now is 03:36.