[kalaikdr]

Hi,

I am trying to implement the CLSVOF scheme in interFoam. as given to Albadawi et al, International Journal of Multiphase Flow, 2013.

In his paper, it is given that the redistancing of the level set function is solved by using the equation

d(phi)/dt' = sign(phi)*(1-|grad(phi)|)

where t' is an artificial time step chosen as delta t' = 0.1*delta x

when I tried it on interFoam as

solve
(
fvm::ddt(lsPhi)
+(sign(lsPhi) * mag(fvc::grad(lsPhi)))
-sign(lsPhi)
);

Foam compiles, but crashes on running due to dimension mismatch.

Any idea how to set it right? Does t' have the units of distance, as it is based on the mesh size?

Also, how do I solve this in a loop of artificial time until it converges to the steady state solution?

[feddy]

Hi
Try this

1- you define the field d and d0 (d is your lsPhi)

volScalarField d
(
IOobject
(
"d",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);

volScalarField d0
(
IOobject
(
"d0",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("d0",dimensionSet(0,0,0,0,0,0,0) ,scalar(0)),
d.boundaryField().types()
);

2- d have no dimension, for do not have one error of dimension you define:

dimensionedScalar one = dimensionedScalar("one", dimensionSet(0,1,0,0,0,0,0), 1.0);

3- you define the quantity present in the article:

const surfaceScalarField deltaCoeff = d.mesh().deltaCoeffs();
const scalar DX = (1.0/min(deltaCoeff).value());

scalar Dtau = 0.1*DX;
scalar Gamma = 0.75*DX;
scalar EPS = 0.1*DX;
scalar NITER = EPS/Dtau;

4- the pseudo iteration

d == Gamma*(2*alpha1-1.0);
d0 == d;

for (int loop = 0; loop < NITER; loop++)
{
// calculate distance function field

fvScalarMatrix dEqn
(
fvm::Sp(scalar(1),d) + Dtau*sign(d0)*mag(one*fvc::grad(d)) == d + Dtau*sign(d0)
);
dEqn.solve();
}

[P Sharma]

Hii everyone,
I am new in openFoam. I want to know how to use CLSVOF model in openFoam (or in foam extend) . Or I have to implement clsvof model separately. Is their any chance to install that model direct and i can use.

[MajidHagh]

Albadawi claimed that he gets less spurious currents with s-clsvof. But I've tried every thing he mentioned and double checked my solver times. I really doubt about his results since in his next paper he didn't use s-clsvof whereas he showed great result in this paper.
Anyone tried s-clsvof to get lower spurious current around the interface?

[mheinz]

I am going to try the same thing within the next one or two weeks. If I find something, I will let you know.

[Saideep]

Hi guys;

i have a doubt regarding VOF's alpha. Interfoam's alpha calculated in some cases (at lower Capillary number flows, surface tension dominated flows) is highly smeared out. In such cases how can you link alpha and phi? There is huge error coming from alpha term everything else seems to work perfectly{like the calculation of surface tension force, capillary pressure etc}.

Any advice over this issue,
Saideep

[Saideep]

Hi Foamers,

Any development with this solver? The surface tension forces seen to be better than interfoam calculations but the velocities at the interface are shooting up.

Meanwhile I have been trying to understand an other parallel solver named "VOSET".

http://ac.els-cdn.com/S0045793015003...ea63ef66bfe6ce

It looks much complicated and I think implemented in Fluent and not in OpenFOAM.

Any share of updates;
Saideep

[mheinz]

Hello,

my experience with the S-CLSVOF solver so far is, that the solution is very sensitive to the value specified for the mesh step size (DX in the above post of feddy).

In my case, I implemented the energy equation in interFoam. When I utilize the Heaviside function to update all the fluid properties and the face fluxes, the temperature gets unbounded after some time steps.

Just physically, the results seem to be ok when I stick to updating the fluid properties and the face fluxes with alpha. However, I still need to compare some results with relevant test cases to see quantitative relations. Furthermore I will implement the direct specification of the mesh step size for the relevant regions instead of setting one fix value.

Regarding the other method you just mentioned: I did not take a closer look to the paper but the coupling seems to be done fully geometrically. As far as I can say, the geometrical reconstruction of the interface is always the more accurate method but comes with a substantial coding effort. Therefore it would depend on your resources and coding skills whether you want to implement something like this.

Regards,
Michael

[Saideep]

Hi,

but how is your velocity behavior? Do you observe any shoot ups at the interface surface or does it stabilize to an average value over time?

By the way what Capillary number are you dealing with?

Saideep

[mheinz]

Hiho,

so far I did not experience odd velocities at the interface.

My capillary number goes up to approximately 0.1, rising along with the runtime.

Regards,
Michael

[Saideep]

Hi,

sorry to get back late.

I am dealing with extremely low Ca values around e-5 or so.

Did you try your code over a static droplet? If yes how does it behave?

In general VOF is reliable for higher Ca flows but for very low Ca flows there seems to be a problem with random velocities at the interface of fluids.

Saideep

[mheinz]

Hi.

I did some further reading on the influence of spurious currents since I am quite new to OpenFOAM and VOF methods in general. My Ca number of 0.1 seems to be indeed quite high compared to cases where spurious currents are strongly recognized.

I did not simulate a static droplet yet. However, I was planning to do this as one validation case of my final model. When I have any results I can send you the details.

Regarding your spiking velocities at the interface - As far as I understand they are a result of the inaccurate curvature calculation from the alpha1 field. How strong did you couple the VOF and LS methods?
Do you just calculate the distance field of the LS method to recalculate the unit normal vectors/ curvature in order to update the fluxes and calculate the surface tension force in UEqn.H and pEqn.H?

You could also consider to update the fluxes and fluid properties with the Heaviside function and use the newly calculated fluxes in UEqn.H. I was also thinking to use the unit normal vector calculated from the LS field for the calculation of phir in the original alphaEqn.H. Considering the last step I am not sure if this will make the initialization of the LS field somewhat corrupted.

Any thoughts?

Regards,
Michael

[P Sharma]

Hi,
anyone knows what is time scale for spurious current in both when we use adaptive grid or uniform grid.

[Saideep]

Hi,

Well, I just followed Albadawi's paper and tried both cases of linking my fluid properties with "Heaviside function" and using VOF's indicator function directly to be applied to compute the fluid properties. Though I found not much of difference.

@ Pravven:
I guess there is no time frame to see the spurious currents evolve. Initially during relaxation step, we see velocities but over time they should decompose, in general they dont decompose but generate vorticities due to force imbalance.

Saideep

[mheinz]

Hi.

I think that is not a bad thing. There should not be big differences. How exactly do you update your fluxes? I am still having problems when I do one recalculation of the alphaEqn with the MULES solver.

Regarding my problems I am thinking about introducing a smoothed alpha field for the calculation of interface properties and just calculate the distance field without any coupling. That is only because I want to use the distance in phase change calculations.

Maybe this would be another approach for your case too?

Regards,
Michael

[Saideep]

Hi Michael,

What do you mean by recalculation of alpha? You are calculation alpha advection only once at the start of the piso/ pimple loop.

I have tried some filtering of alpha and interface curvature etc and seems to be fine for some cases.

I am relatively new to level set method and well i was having the similar idea as you. I don't understand why do we need to run the loop to re-distance the LS function.

If I specify the initialization of the LS function as (2*alpha1 -1) it is already on the interface(alpha1 = 0.5).

From other papers where coupling of Ls with VOF is done, they advect the LS function at the start of the loop so it is essential to re-distance the function. But here we aren't advecting the LS function.

Saideep

[mheinz]

Hi.

Albadawi suggests the following steps, among others:

[...]
5. Re-initialize the LS function using Eq. (14) in order to obtain the
signed distance function and correct the interface at the bound-
aries. Then, calculate the new values of the Heaviside function,
the Dirac functions, and the interface curvature.

6. Update the fluid physical properties and the fluxes using the
volume fraction function a (Eqs. (7) and (8)).

7. Solve the Navier Stokes system of equations for velocity and
pressure using the Pressure Implicit with Splitting of Operators
(PISO) (Issa, 1986). Details of its implementation in OpenFOAM
are provided in the Appendix. In the present study three pres-
sure correction steps were used and ensured that the continuity
residual remained always below 1e-6.

After the calculation of the distance field (#5) he suggests updating the fluxes and fluid properties (#6). At this point I calculate the following:

{
    word alphaScheme("div(phi,alpha)");
    word alpharScheme("div(phirb,alpha)");

    surfaceScalarField phic(mag(phi/mesh.magSf()));
    phic = min(interface.cAlpha()*phic, max(phic));

    #include "calcK.H"
    surfaceScalarField phir(phic*nHatf);
    //surfaceScalarField phir(phic*interface.nHatf());

    surfaceScalarField phiAlpha
    (
        fvc::flux
        (
            phi,
            alpha1,
            alphaScheme
        )
      + fvc::flux
        (
            -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
            alpha1,
            alpharScheme
        )
    );

    MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);

    rho == alpha1*rho1 + (scalar(1) - alpha1)*rho2;

    volScalarField& nu = const_cast<volScalarField&>(twoPhaseProperties.nu()());
    nu == alpha1*nu1 + (scalar(1)-alpha1)*nu2;

    rhoCp == alpha1*rho1*cp1 + (scalar(1)-alpha1)*rho2*cp2;

    rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;

    rhoPhiCp = phiAlpha*(rho1*cp1 - rho2*cp2) + phi*rho2*cp2;
}

This corresponds to one iteration of alphaEqn with the new unit normal vector nHatf, calculated in "calcK.H". At this point one could also use the heaviside function to update the fluxes and fluid properties. For the fluxes one needs to simply change "alpha1" and "phiAlpha" to "heaviside" and "phiHeaviside". Since I am not getting any reasonable results, I think this step is wrong. Maybe it would help to do more iterations, e.g. by adding a alphaEqnSubCycle-like routine. However, since I think the problem in my case originates in the calculation of the interface properties I will not take this approach further for now.

The redistancing is still a necessary step. By initializing of the distance field you just take a "good guess" with the alpha field as a basis for the position of the interface.

psi0 == (2*alpha1 - 1)*epsilon;
psi == psi0;

There you have +epsilon in the liquid phase and -epsilon in the vapor phase. To get a signed distance field, you need to do a correction:

for (int corr=0; corr<int(epsilon.value()/deltaTau); corr++)
{
    psi = psi + sign(psi0)*(scalar(1) - mag(fvc::grad(psi)*dimM))*deltaTau;

    psi.correctBoundaryConditions();
}

Since we use an artificial time step, we do not need to solve a transport equation for the distance field, hence the "simple" coupling.

What kind of smoothing did you use in your case? I am thinking to implement the laplacian filter Ubbink mentions in his work (Onno Ubbink, 1997, Numerical prediction of two fluid systems with sharp interfaces, Eq. 3.51).

Regards,
Michael

[vigneshTG]

Guys,

I tried S-CLSVOFFoam version by Takuya Yamamoto . He has modified the original version by A. Albadawi (i think) where he solves the volume fraction equation with interface normal calculated from the level set function which by default updates the flux using the level set. It appears to be a little better, but as Saideep pointed there are some odd velocities around the interface which seems to increase with mesh refinement. Nevertheless, the accuracy improves as the mesh is refined. (test case: rise of meniscus close to a wall)

He has also another version where he uses different relations for calculating the delta and Heaveside function

[abdoli]

Dear MajidHagh
do you try to regenerate numerical simulation data of bubble formation by air injection as is mentioned by Albadawi et al. (2013).
I want to do that but my results are different. If yes, would you mind send me your setting files (0,constant,system).

thanks a lot
abdoli

[MajidHagh]

I've worked on this method for a year and I faced problems as you did.
Recently I improved the method and I'm working on a paper which describes the correct way of implementing s-CLSVOF. I've simulated couple of two-phase cases to evaluate my solver and I got good result.
I would release my solver online when I publish the paper. Hopefully it would be soon.
Majid