[Andrea_85]

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

[akidess]

hg clone https://code.google.com/r/robertocastillalopez-interfoamssf-210/

[lindstroem]

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:

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

[akidess]

The output of foamCalcEx looks like this:

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

[lindstroem]

ah ok, sorry, now i got how your script was supposed to run

Thanks!

[Andrea_85]

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

[akidess]

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

[Andrea_85]

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

[akidess]

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

[Andrea_85]

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]

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]

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

[akidess]

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

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

[Andrea_85]

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

[akidess]

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)!

[Andrea_85]

good job! you are right.

andrea

[Andrea_85]

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

[akidess]

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?

[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

[akidess]

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

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