CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

VOF interFoam - Strange behaviour in a micro T-junction

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

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   November 17, 2015, 06:27
Default VOF interFoam - Strange behaviour in a micro T-junction
  #1
Member
 
Paolo Capobianchi
Join Date: Sep 2014
Posts: 35
Rep Power: 9
pablitobass is on a distinguished road
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
Attached Images
File Type: jpg drop01.jpg (18.8 KB, 144 views)
File Type: jpg drop02.jpg (19.7 KB, 113 views)
File Type: jpg drop03.jpg (20.5 KB, 102 views)
File Type: jpg drop04.jpg (21.2 KB, 100 views)
File Type: jpg drop05.jpg (21.3 KB, 109 views)
pablitobass is offline   Reply With Quote

Old   November 17, 2015, 06:59
Default
  #2
Senior Member
 
akidess's Avatar
 
Anton Kidess
Join Date: May 2009
Location: Germany
Posts: 1,377
Rep Power: 28
akidess will become famous soon enough
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.
__________________
*On twitter @akidTwit
*Spend as much time formulating your questions as you expect people to spend on their answer.
akidess is offline   Reply With Quote

Old   November 17, 2015, 07:22
Default
  #3
Member
 
Paolo Capobianchi
Join Date: Sep 2014
Posts: 35
Rep Power: 9
pablitobass is on a distinguished road
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:

Code:
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
pablitobass is offline   Reply With Quote

Old   November 17, 2015, 07:32
Default
  #4
Senior Member
 
akidess's Avatar
 
Anton Kidess
Join Date: May 2009
Location: Germany
Posts: 1,377
Rep Power: 28
akidess will become famous soon enough
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.
__________________
*On twitter @akidTwit
*Spend as much time formulating your questions as you expect people to spend on their answer.
akidess is offline   Reply With Quote

Old   November 17, 2015, 07:57
Default
  #5
Member
 
Paolo Capobianchi
Join Date: Sep 2014
Posts: 35
Rep Power: 9
pablitobass is on a distinguished road
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?

Last edited by pablitobass; November 17, 2015 at 10:30.
pablitobass is offline   Reply With Quote

Old   November 18, 2015, 06:01
Default
  #6
Member
 
Paolo Capobianchi
Join Date: Sep 2014
Posts: 35
Rep Power: 9
pablitobass is on a distinguished road
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
pablitobass is offline   Reply With Quote

Old   November 19, 2015, 13:44
Default
  #7
Senior Member
 
akidess's Avatar
 
Anton Kidess
Join Date: May 2009
Location: Germany
Posts: 1,377
Rep Power: 28
akidess will become famous soon enough
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.
__________________
*On twitter @akidTwit
*Spend as much time formulating your questions as you expect people to spend on their answer.
akidess is offline   Reply With Quote

Old   November 20, 2015, 04:59
Default
  #8
Member
 
Paolo Capobianchi
Join Date: Sep 2014
Posts: 35
Rep Power: 9
pablitobass is on a distinguished road
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
pablitobass is offline   Reply With Quote

Old   November 20, 2015, 05:32
Default
  #9
Senior Member
 
akidess's Avatar
 
Anton Kidess
Join Date: May 2009
Location: Germany
Posts: 1,377
Rep Power: 28
akidess will become famous soon enough
Yeah, I never run a simulation without momentum predictor. Which does lead to the question though, how are your pressure residuals behaving?
__________________
*On twitter @akidTwit
*Spend as much time formulating your questions as you expect people to spend on their answer.
akidess is offline   Reply With Quote

Old   November 20, 2015, 06:13
Default residuals
  #10
Member
 
Paolo Capobianchi
Join Date: Sep 2014
Posts: 35
Rep Power: 9
pablitobass is on a distinguished road
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
Attached Images
File Type: jpg residuals_moment_predict_off.jpg (57.0 KB, 84 views)
File Type: jpg residuals_moment_predict_on.jpg (79.3 KB, 69 views)

Last edited by pablitobass; November 20, 2015 at 10:08.
pablitobass is offline   Reply With Quote

Old   November 20, 2015, 06:51
Default
  #11
Senior Member
 
akidess's Avatar
 
Anton Kidess
Join Date: May 2009
Location: Germany
Posts: 1,377
Rep Power: 28
akidess will become famous soon enough
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.
__________________
*On twitter @akidTwit
*Spend as much time formulating your questions as you expect people to spend on their answer.
akidess is offline   Reply With Quote

Old   November 20, 2015, 07:26
Default
  #12
Member
 
Paolo Capobianchi
Join Date: Sep 2014
Posts: 35
Rep Power: 9
pablitobass is on a distinguished road
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 is offline   Reply With Quote

Old   November 20, 2015, 07:45
Default
  #13
Member
 
Paolo Capobianchi
Join Date: Sep 2014
Posts: 35
Rep Power: 9
pablitobass is on a distinguished road
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 is offline   Reply With Quote

Old   November 21, 2015, 15:18
Default
  #14
Member
 
Paolo Capobianchi
Join Date: Sep 2014
Posts: 35
Rep Power: 9
pablitobass is on a distinguished road
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
Attached Images
File Type: jpg residuals.jpg (63.7 KB, 33 views)
File Type: jpg drp_detachment.jpg (22.2 KB, 71 views)
pablitobass is offline   Reply With Quote

Old   November 22, 2015, 10:33
Default
  #15
Senior Member
 
akidess's Avatar
 
Anton Kidess
Join Date: May 2009
Location: Germany
Posts: 1,377
Rep Power: 28
akidess will become famous soon enough
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.
__________________
*On twitter @akidTwit
*Spend as much time formulating your questions as you expect people to spend on their answer.
akidess is offline   Reply With Quote

Old   November 22, 2015, 11:30
Default
  #16
Member
 
Paolo Capobianchi
Join Date: Sep 2014
Posts: 35
Rep Power: 9
pablitobass is on a distinguished road
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:

Code:
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

Code:
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
pablitobass is offline   Reply With Quote

Old   November 24, 2015, 02:06
Default
  #17
Senior Member
 
akidess's Avatar
 
Anton Kidess
Join Date: May 2009
Location: Germany
Posts: 1,377
Rep Power: 28
akidess will become famous soon enough
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.
__________________
*On twitter @akidTwit
*Spend as much time formulating your questions as you expect people to spend on their answer.
akidess is offline   Reply With Quote

Old   November 24, 2015, 05:38
Default
  #18
Member
 
Paolo Capobianchi
Join Date: Sep 2014
Posts: 35
Rep Power: 9
pablitobass is on a distinguished road
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
pablitobass is offline   Reply With Quote

Old   December 7, 2015, 07:09
Default
  #19
Senior Member
 
akidess's Avatar
 
Anton Kidess
Join Date: May 2009
Location: Germany
Posts: 1,377
Rep Power: 28
akidess will become famous soon enough
Any luck?
__________________
*On twitter @akidTwit
*Spend as much time formulating your questions as you expect people to spend on their answer.
akidess is offline   Reply With Quote

Old   December 7, 2015, 16:30
Default
  #20
Member
 
Paolo Capobianchi
Join Date: Sep 2014
Posts: 35
Rep Power: 9
pablitobass is on a distinguished road
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
pablitobass is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
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 Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Strange behaviour when using compressibleInterFoam with constantAlphaContactAngle TobM OpenFOAM Running, Solving & CFD 2 May 11, 2016 06:34
Strange temperature behaviour with interFoam tayo OpenFOAM 9 July 11, 2014 07:28
strange pressure behaviour with symmetricPlane boudary condition - interFoam duongquaphim OpenFOAM Running, Solving & CFD 10 August 20, 2013 14:00
Strange boundary behaviour using interFoam Andrea_85 OpenFOAM 11 January 22, 2013 15:09
Strange behaviour because of contact angle (interfoam) Kim123 OpenFOAM Running, Solving & CFD 0 January 12, 2011 10:16


All times are GMT -4. The time now is 21:56.