[pablitobass]

Hi All,

I am trying to simulate a flow of two immiscible liquids in a micro T-juntion (characteristic length 100 micron) using the interFoam solver and the results I get are probably wrong.

In the flow regime I am studying, I should see the detachment of droplets in the main channel and the point of detachment must be fixed in space. However, in my case the drops detachment point "walks" trough the channel (the attachments show only 5 droplets, but I have run the simulation for much longer)

I guess that it is due to the spurious currents, and I have tried many different set-up to improve the results. I did the following things:

- tried different level of refinement of the mesh (coarser meshes should be less problematic for parasite currents);
- I have used a "smoothed" volume fraction: I interpolate alpha in the cell faces and then I calculate the average over the entire surface of the cell.I tried to apply 1, 2 and 3 level of smoothing (the "recommended" is 2 times);
- I have used a more accurate scheme for the calculation of the gradient of alpha (extendedLeastSquares 0.0);
- I have used "cubic" as default in the interpolationSchemes, although I am not sure that it can helps;
- I have used a maximum of 10 alpha subcycles and 5 alpha correctors in order to have a "reasonable" time step and accuracy for the calculation of the volume fraction;
- the maxCo for the momentum equation has been set to 0.5 and 1.

To be honest, I do not know if I can do something else/more to improve the results. If someone has any suggestion, would be really great!

Thanks very much in advantage!

Regards,
Paolo

[akidess]

Are your simulations 2D? That might not reflect the physics involved. Other than that, maybe your boundary conditions are wrong (main channel should be hydrophobic?), or you need to further refine your mesh. Spurious currents is *not* the first thing that comes to mind.

[pablitobass]

Hi Anton,

Thank you very much for your reply!

The case is 3D, I forgot to mention. The images are slices. About the boundary conditions, I have set hydrophobic conditions for the dispersed phase. This is my alpha1 file settings:

internalField   uniform 0;

boundaryField
{
    wall
    {
        type             constantAlphaContactAngle;
        limit              gradient;
        theta0          180;
        value            uniform 0;
    }

    inlet_D_wall
    {
        type             constantAlphaContactAngle;
        limit             gradient;
        theta0          0;
        value           uniform 1;
    }

    inlet_C
    {
        type            fixedValue;
        value          uniform 0;
    }

    inlet_D
    {
        type             fixedValue;
        value           uniform 1;
    }

    outlet
    {
        type            zeroGradient;
    }
}

About the mesh, I tried a quite refined one, but the problem is still there.

Cheers,
Paolo

[akidess]

So the red phase should never touch the wall, right? Well it does seem to do so in figure 5 (by the way, don't look at Paraview smoothed fields when trying to fix problems). Perhaps fixedValue 0 works better.

[pablitobass]

Yes Anton, the red fluid is the dispersed phase that should not touch the wall. However, as you have noticed, as the dispersed phase advance along the channel, the phase seems to wet the wall.Maybe it is due to Paraview?

So, you suggest something like this for the dispersed phase:

type constantAlphaContactAngle;
theta0 fixedValue 0;
value uniform 1;

Am I wrong?

[pablitobass]

Just few additional informations. I did simulations for lower viscosities for the continuous phase (the blue fluid) and it works fine.

I do not really know why for higher viscosities it does not work. In principle, if we assume that the problem is due to the spurious currents, lower capillary numbers, i.e. lower viscosities of the continuous phases, should be less problematic.

To be honest, I do not have any clue. Maybe it is a problem of boundary conditions? But if it so, why the lower viscosity cases works fine (I have compared the result with other taken from a paper)?

These simulations are fundamental for my PhD, and I really have to find out if the problem is due to my set up or not.

Thanks a lot in advantage!

Regards,
Paolo

[akidess]

AFAIK the value of the viscosity is less important than the viscosity ratio. If you increase the viscosity ratios, you will run into more problems.

In your case, I wouldn't use a contact angle boundary condition. For fully wetting or non-wetting conditions, I'd just use a regular fixedValue boundary condition with value 0 or 1. For 90 degree contact angle I'd use zeroGradient.

By the way, if your domains are always this simple you might also want to look at more sophisticated codes, such as Gerris or Paris.

[pablitobass]

Hi Anton,

Thank you very much for you reply.

I agree with you that the viscosity ratio has more impact on the solution instead of the value for a single phase. When I said I have increased the viscosity of the continuous phase I meant implicitly I have increased the viscosity ratio.

In the meanwhile, I had a deep look at the discussion about interFoam and I found something interesting. It seems that other users had problem with large viscosity ratios and a possible remedy would be to switch the momentum predictor on. I tried to do so, and I have noticed an improvement although the detachment still moves, but "slower".

I have just run a new simulation whit the boundary condition you suggested and the same set up as before. Tomorrow I will have an idea if it will improve the solution. I will let you know.

I will have a look at the solvers you suggested, although I am afraid that at this point it is too late to change for me. I have to work also with non-Newtonian fluids, and it seems that they do not support that kind of models.

Thanks again!

Regards,
Paolo

[akidess]

Yeah, I never run a simulation without momentum predictor. Which does lead to the question though, how are your pressure residuals behaving?

[pablitobass]

Hi Anton,

I have attached two pictures with the residuals, with (second image) and without momentum predictor. I might be wrong, but in both cases they don't seem very good...

Regards,
Paolo

[akidess]

I agree, your residuals seem too large. What's your time step and Courant number (max)? Compare it with the stability limit from Deshpande et al, 2012 ( doi:10.1088/1749-4699/5/1/014016 ). Decrease your time step and increase your pressure correctors.

[pablitobass]

Hi again,

my max Co is about 0.5 for both momentum and volume fraction equations in the current simulation. The time step is about 4e-07 which is lower than the value I have calculated previously accordingly to the paper you suggested me and the paper of Galusinsky and Vigneaux.

Now I have increased the number of correctors and decreased the max Co. Will see what will happen.

Thanks for the suggestions.

Regards,
Paolo

[pablitobass]

Just a correction to the previous post. The value I had calculated previously was based on another degree of refinement of the mesh. I have just recalculated it and my time step should be less than 1.5e-07. I have changed the settings accordingly.

Regards,
Paolo

[pablitobass]

Hi Anton,

I tried to do what you suggested. I have changed the boundary conditions for aplha, decreased the time step and increased the number of correctors.
Now the dispersed phase does not touch the wall at all and the residuals look better, although are still too high probably.

Unfortunately, the problem for the detachment point is still there. I have no clue on how to solve it, if we admit that can be solved...

Regards,
Paolo

[akidess]

What's the value of cAlpha you are using? You have some spurious satellite droplets behind the bigger ones. Also post your fvSchemes. Apart from that, I'm running out of ideas too unfortunately.

[pablitobass]

Hi Anton,

Thank you for your reply!

I have used cAlpha = 0.75 because previously I thought the problem was due the spurious currents, so a bit of "interface smearing" maybe could help. Now I am running two jobs with cAlpha = 1. I have also increased the number of pressure correctors and further reduced the time step. Now the residuals are much better, but unfortunately the speed of calculation has reduced enormously.

This is the fvSchemes:

ddtSchemes
{
    default         CrankNicolson 1;
}

gradSchemes
{
    default         extendedLeastSquares 0.0;
}

divSchemes
{
    div(rho*phi,U)  Gauss limitedLinearV 0;
    div(phi,alpha)  Gauss vanLeer 01;
    div(phirb,alpha) Gauss interfaceCompression;
    div((muEff*dev(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         cubic; 
}

snGradSchemes
{
    default         corrected;
}

fluxRequired
{
    default         no;
    p_rgh;
    pcorr;
    alpha1;
}

and this is the fvSoultion

solvers
{
    pcorr
    {
        solver          PCG;
        preconditioner
        {
            preconditioner  GAMG;
            tolerance       1e-10;
            relTol          0;
            smoother        GaussSeidel;
            nPreSweeps      0;
            nPostSweeps     2;
            nBottomSweeps   2;
            cacheAgglomeration false;
            nCellsInCoarsestLevel 10;
            agglomerator    faceAreaPair;
            mergeLevels     1;
        }

        tolerance       1e-8;
        relTol          0;
        maxIter         100;
    }

    p_rgh
    {
        solver          GAMG;
        tolerance       1e-08;
        relTol          0.01;
        smoother        GaussSeidel;
        nPreSweeps      0;
        nPostSweeps     2;
        nFinestSweeps   2;
        cacheAgglomeration false;
        nCellsInCoarsestLevel 10;
        agglomerator    faceAreaPair;
        mergeLevels     1;
    }

    p_rghFinal
    {
        solver          PCG;
        preconditioner
        {
            preconditioner  GAMG;
            tolerance       1e-08;
            relTol          0;
            nVcycles        2;
            smoother        GaussSeidel;
            nPreSweeps      0;
            nPostSweeps     2;
            nFinestSweeps   2;
            cacheAgglomeration false;
            nCellsInCoarsestLevel 10;
            agglomerator    faceAreaPair;
            mergeLevels     1;
        }

        tolerance       1e-08;
        relTol          0;
        maxIter         100;
    }

    U
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-08;
        relTol          0;
        nSweeps         1;
    }

    UFinal
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        preconditioner  GAMG;
        tolerance       1e-08;
        relTol          0;
        nSweeps         1;
    }
}

PIMPLE
{
    momentumPredictor yes;
    nCorrectors     15;
    nNonOrthogonalCorrectors 0;
    nAlphaCorr      10;
    nAlphaSubCycles 10;
    cAlpha          0.75;

}

Thanks again for you effort on helping me!

Regards,
Paolo

[akidess]

Since you are debugging, I'd recommend using an Euler time marching scheme and a linear interpolation scheme (ideally test the effect of these settings using two simulations). Then I think you are trying to use the bounded version of vanLeer, which should be vanLeer01 (without a space!). Finally, I'd recommend trying OF-3.0.x, since there has been some work on the multiphase model.

[pablitobass]

Hi Anton,

Thanks for the reply.

I will try the 3.0.x version, maybe works better.

About the schemes, I have tried many different simulation with the Euler time scheme and the linear interpolation as it was by default. I have changed them to see if could solve the issue of the detachment point advancing.

About the vanLeer01, thanks for the tip, I did not noticed that. In practice, I guess I was working with an unbounded scheme. Maybe that's why I had so much smearing of the interface.

Thanks once more for your advices!

Regards,
Paolo

[akidess]

Any luck?

[pablitobass]

Hi Anton!

Thanks for asking. Actually is getting better... Now I am using interDyMFoam and when the simulation reach a steady configuration the pinch-off point does not move anymore. Also the diameter of the droplets is in good agreement with the results from the reference paper. It seems that was a problem of the mesh, however, before I could not run at all a case without the adaptime mesh with the same degree of refinement...

On the other hand, there are still some problems. There is a not negligible smearing of the interface and the droplets have some unrealistic "white spots" due to the inclusion of the continuous phase during the break up. Maybe it is due to the choice of the flux limiter for the advective term (vanLeer01), but I am not sure though.

Furthermore, consider that because of the large increasing of the domain size, I could not set a time step that guarantee the stability condition imposed by the new mesh size at the interface. The pressure residuals look not really well. Maybe it can have an impact as well on the above mentioned problem. I am thinking about a local time stepping...

This is pretty much the sum of what I have done, but if you have any suggestion, would be rellay appreciated

Regards,
Paolo