[Andrea_85]

HI all,

Does force calculation work for multiphase cases?. If i try to add these lines below to my controldict in damBreak example:

functions
(
forceA
{
type forces;
functionObjectLibs ("libforces.so");
patches (leftWall);
rhoInf 1000;
CofR (0 0 0);
outputControl timeStep;
outputInterval 1;
pName p;
UName U;
rhoName rho;
log true;
}


i got the following error:

keyword nu is undefined in dictionary "/home/aferrari/OpenFOAM/OpenFOAM-2.1.0/tutorials/multiphase/interFoam/laminar/damBreak/constant/transportProperties"

file: /home/aferrari/OpenFOAM/OpenFOAM-2.1.0/tutorials/multiphase/interFoam/laminar/damBreak/constant/transportProperties from line 20 to line 62.

From function dictionary::lookupEntry(const word&, bool, bool) const
in file db/dictionary/dictionary.C at line 400.

FOAM exiting


I have already serched the forum and found how to get around the problem, but the question is if force calculation is suitable for multiphase simulations or not. What does it mean for example rhoInf in this case? rho is not constant in space. Even the viscosity is the viscosity of the mixture calculated in twoPhaseMixture.C, not simply mu_1 or mu_2.

Can someone explain this?

regards

andrea

[niaz]

dear andrea
if you look at force calculation you can find it.
it is not suitable for this purpose.
you should develop a code to save nu as mixed of two phase and
convert it to single phase.

[Andrea_85]

Thanks niaz,
I thought it was not suitable for this, just to be sure.
Now what i want to do is try to calculate pressure and viscous forces
separately in both phases (maybe based on alpha). Does it make sense in your opinion?

regards

andrea

[calim_cfd]

Originally Posted by Andrea_85

Thanks niaz,
I thought it was not suitable for this, just to be sure.
Now what i want to do is try to calculate pressure and viscous forces
separately in both phases (maybe based on alpha). Does it make sense in your opinion?

regards

andrea

sorry to interfere but i guess that first you should check what flow pattern your working with for this multiphase flow... then i guess you should be able to decide whether an average value or separated values should suit your application best,

cuz pressure is just the normal forces and viscous the tangential ones

if this does not make sense i'm sry and pls correct me

as far as the application goes i guess it should be straightforward if ur familiar with C++

regards
/calim

[niaz]

Andrea
if you look at base file for force and forcecoeff you can find that both are for
single-phase flow. but in interdyfoam tut, a sample used it for floating object and it has answered well.
i recommend to use it for a standard test case to be sure about it.

[Andrea_85]

Hi Niaz and Calim,
thanks for your suggestions. This is what i've done:

in forces.C where devRhoReff() is calculated i added:

dimensionedScalar nu1(transportProperties.subDict("phase1").lookup("nu"));
       dimensionedScalar nu2(transportProperties.subDict("phase2").lookup("nu"));
       dimensionedScalar rho1(transportProperties.subDict("phase1").lookup("rho"));
       dimensionedScalar rho2(transportProperties.subDict("phase2").lookup("rho"));

const volScalarField limitedAlpha1
     (
         min(max(alpha1, scalar(0)), scalar(1))
     );
        
            
     const volScalarField mu
         (
     
             limitedAlpha1*rho1*nu1
            + (scalar(1) - limitedAlpha1)*rho2*nu2
         );

return mu*dev(twoSymm(fvc::grad(U)));

and in calcForcesMoment()

const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);
const volScalarField& p = obr_.lookupObject<volScalarField>(pName_);
const volScalarField& alpha1 = obr_.lookupObject<volScalarField>("alpha1");

const fvMesh& mesh = U.mesh();
     
        const surfaceVectorField::GeometricBoundaryField& Sfb =
            mesh.Sf().boundaryField();
            
            
        tmp<volSymmTensorField> tdevRhoReff = devRhoReff();
     const volSymmTensorField::GeometricBoundaryField& devRhoReffb
            = tdevRhoReff().boundaryField();
                        
        const volScalarField limitedAlpha1
        (
           min(max(alpha1, scalar(0)), scalar(1))
        );

forAllConstIter(labelHashSet, patchSet_, iter)
        {   
            label patchi = iter.key();

            
vectorField pf1(Sfb[patchi]*(p.boundaryField([patchi]*limitedAlpha1.boundaryField()[patchi]));
vectorField pf2(Sfb[patchi]*(p.boundaryField()[patchi]*(scalar(1)-limitedAlpha1.boundaryField()[patchi])));


vectorField vf1((Sfb[patchi]*limitedAlpha1.boundaryField()[patchi]) & devRhoReffb[patchi]);
vectorField vf2((Sfb[patchi]*(scalar(1)-limitedAlpha1.boundaryField()[patchi])) & devRhoReffb[patchi]);



fmPressurePhase1 += sum(pf1);
fmPressurePhase2 += sum(pf2);


fmViscousPhase1 += sum(vf1);
fmViscousPhase2 += sum(vf2);

and some other changes to print them out. Don't know if it makes sense, i need to test it. I think a standard test case could be poiseuille flow between two parallel plates with 2 fluid and a meniscus between them. When there is only phase 1 (at the beginning) or 2(at the end of sim) i should find analytical solution for viscous dissipation.
What do you think?

best

andrea

[daveatstyacht]

Andrea,
I think you may be double accounting in the case of turbulence models because the effects of variations of nu are captured in the devReff of the individual turbulence models (see kOmega.C for example). The variations in density should be accounted for by the volScalarField rho() if the rhoName is set to the name of rho for that solver (typically just "rho"). Setting rhoName to rhoInf will result in the use of a uniform density volScalarField. Your approach may be necessary for laminar cases, and/or where a breakdown of forces in each phase is desired. Your approach also has the effect of treating faces as either completely one phase or the other which may under-predict forces for flows with large mixture regions.

Dave

[Andrea_85]

Hi Dave,
you are right but i'm working with multiphase flow and my cases are always laminar (0.01<Re<1), for turbulence cases the approach has to be corrected. What i have done is basically create a new "forcesMultiPhase" dir in which i have only what i need.
Why you are saying that i'm treating faces "as either completely one phase or the other"?. the face area vectors are multiplied by alpha so if alpha=1 the face belongs to phase 1, if alpha=0 the face belongs to phase2 and at the intereface the face belongs to both depending on alpha. Can you please explain this?

best

andrea

[daveatstyacht]

Andrea,
I was mistaken reading the setting of limitedAlpha1 and the ordering of the min and max used (I thought it was forcing intermediary alphas to be 0). Since you are working with exclusively laminar flow your approach makes more sense now. I thought of this issue myself, but since I work with fully turbulent flow I concluded it wasn't necessary for my own flow type (ships).

Dave

[Andrea_85]

Hi again,
i would like to share my test case. I finally decided to test force calculation on this case: a slug in a vertical capillary tube with fixed contact angle that moves under gravity. The falling velocity is constant (after an initial acceleration) and the forces acting on the slug should be gravity and viscous forces only (pressure and capillary forces should be zero for symmetry). I tested 2 cases: one in which the gravity is set to (0 -0.1 0) and the other one in which gravity is set to its normal value (0 -9,81 0). the case is 2D. The calculated viscous forces balances the gravity force only in the first case, when the gravity is lower. In case of higher gravity the curves do not overlap (maybe there are other effects due to time derivative..). The calculated viscous force on phase 1 is oscillating a lot (also the total viscous force is oscillating, so i think is not a problem of my implementation) and i do not know exactly why. In the plot attached the black curve is what i called in the previous post "fmViscousPhase1" and the blue line is the gravity force acting on the slug.
Any comments will be appreciated.

best regards

andrea

[calim_cfd]

maybe the higher gravity is changing your flow pattern and the slug pattern is no longer physically possible/expected.. try increasing your volumetric flux, but you might have to deal with turbulence

[Andrea_85]

Hi mauricio

that's what I thought. The Re is around 1600 for the higher gravity case and maybe this could cause the shift. and what's your opinion about the oscillation? The velocity field is more or less constant through the domain except in the region of the interface and i noticed that using a lower courant number i can reduce the magnitude of these oscillations so it might be related to spurios currents.

best

andrea

[calim_cfd]

about parasitic currents.. i'm not that familiar with these cuz i havent rly tried to simulate these flow patterns in cfd.. there dimensionless numbers to check..like capilarity, laplace number..weber's....and i've only had some work with 1d models. but maybe this (http://www.google.com.br/url?sa=t&rct=j&q=spurious+currents+cfd&source=web& cd=4&ved=0CEwQFjAD&url=http%3A%2F%2Fpages.csam.mon tclair.edu%2F~yecko%2Ficodes%2FScardoZaleski_ARFM. pdf&ei=qGRnT4fBK8iN4gS-j_mMCA&usg=AFQjCNGL9JUHaInCbgwZEBAHplLYEiWiMQ) pdf can help you out since you're more into the subject than me... it has some interesting stuff, at least for me.

about the oscillations... it's kinda hard to tell cuz i dont rly know how you're modeling the slug.. i mean..there are a few flow conditions to make it stable. Like for instance the distance from the inlet. you know that a slug only develops at a certain distance from the inlet..say 8 or 16D..so maybe you need a longer/higher tube to get the slug.. plus idk how realistic is the shape of your slug...morever your running a 2d.. i take it's a 2d in cartesian coordinates? anyway i kinda need more time, and this is sth we all lack these days, but if sth comes to mind i'll post it here

[Andrea_85]

HI again mauricio,

I'm really familiar with that paper but i think, at moment, spurios current in interFoam is something you have to "live with". Anyway i uploaded the case if you want to give a try. I checked all forces and found they are all oscillating. Normally i work with bigger domain and the effect of spurios current are only in the time step of sim, they usually do not affect my results. So i do not want to spend much time on this, i was just curios about the reason behind these oscillations.


ps. if you want to try this case just type blockMesh, setField, interFoam, as usual. The sim is a little bit slow, maybe it is better to switch to GAMG.

best
andrea

[calim_cfd]

hi andrea!

thx 4 sharing!
i'll give it a try and let you know if sth comes up

[calim_cfd]

hey andrea i gave it a try but i remembered i dont have your code for force calc.

but i did change sth in your case i did this to your p_rgh file:

 inlet
    {
        type            totalPressure;//fixedValue;
        p0            0;
        rho            rho 1;
        gamma        1;
        value           uniform 0;
    }
    outlet
    {
        type            totalPressure;//fixedValue;
        p0            0;
        rho            rho 1;
        gamma        1;
        value           uniform 0;
    }

cuz as the slug goes down air is gotta go in.. so i set it to total pressure so it corrects the velocity. it's an equivalent bc i guess, and the results seem equal to your orig. settings too

plus i've set gravity to 9.81/m/s2

the slugs seems go down quite smoothly ..so far im unable to tell where to oscillations come from

EDIT: have you tried refining your mesh? The oscilation is very consistent, the case is stable and the slug has a ,say, expected behavior. Then pressure and velocity variations along the wall might be varying too fast for your force calculation due to the higher gravity force.. if you havent tried refining your mesh i rly would like to see the force calculation for a finer one. if the oscillations decrease then it could have been it..idk

then i guess that running it with a lower Co, say, .2 might decrease the oscillations

EDIT2:i've checked the wallstress with oF's appl, assuming air, and i see ripples on the walls that also seem to slide down the walls, but in a oscillating fashion.. even with a finer mesh..that might be the source of the oscillations

here's a pic

1) assuming it's a physical behavior, maybe the system's rly become too dynamic and it's damping is not enough to make things stable..

2) if it's a numeric error i can't tell where it is coming from

[Andrea_85]

Hi mauricio,
thanks for your interest in my post. Yes, i tried to refine the mesh but with no succes in reduce the oscillation. Even use totalPressure instead of fixedPressure bc does not seem to solve the problem. Running the simulation with gravity sets to 9.81 m/s^2 produces a much more smooth curve for viscous force: the oscillation are still there but some order of magnitude lower than the force itself. My guess is that with small velocities the calculated force has the same order of magnitude of the oscillations and you can not filter them out.

Looking at your picture one thingh that comes in my mind is that this behaviour might be related to the formation of a thin flim along the wall, which is an expected/physical behavior. I attacched the "rescaled" color function in which the film is well visible.

best

andrea

[calim_cfd]

Originally Posted by Andrea_85

....My guess is that with small velocities the calculated force has the same order of magnitude of the oscillations and you can not filter them out.
...
Looking at your picture one thingh that comes in my mind is that this behaviour might be related to the formation of a thin flim along the wall, which is an expected/physical behavior. I attacched the "rescaled" color function in which the film is well visible.

yes i see the film too..!

as for the force maybe you cannot scape the scale issue when using a higher gravity.

then i guess that in the end there's nothing wrong with the oscillations? it's just physics (?)

i was interested in this case cuz as i told b4 i recently been through this topic and it's always nice to learn sth new!

[winden]

Dear all.

I don't know if this helps but I was also confused that the solver did not recognize nu,rho or transportModel from my transportProperties. I'm using v.1.7.1 and, looking at the code, this does not seem to matter. I simply defined those three keywords outside of the phase-dictionaries in transportProperties to make it run without crashing and it works fine. To verify that it is actually recognizing the difference in rho and nu, I made the force-library print out the values on the entire patch for plotting instead of just summing them up. See the attached image, same results for pressure force.

//Björn

[niaz]

Dear Andrea
can you send your library for me?

Originally Posted by Andrea_85

Hi Niaz and Calim,
thanks for your suggestions. This is what i've done:

in forces.C where devRhoReff() is calculated i added:

dimensionedScalar nu1(transportProperties.subDict("phase1").lookup("nu"));
       dimensionedScalar nu2(transportProperties.subDict("phase2").lookup("nu"));
       dimensionedScalar rho1(transportProperties.subDict("phase1").lookup("rho"));
       dimensionedScalar rho2(transportProperties.subDict("phase2").lookup("rho"));

const volScalarField limitedAlpha1
     (
         min(max(alpha1, scalar(0)), scalar(1))
     );
        
            
     const volScalarField mu
         (
     
             limitedAlpha1*rho1*nu1
            + (scalar(1) - limitedAlpha1)*rho2*nu2
         );

return mu*dev(twoSymm(fvc::grad(U)));

and in calcForcesMoment()

const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);
const volScalarField& p = obr_.lookupObject<volScalarField>(pName_);
const volScalarField& alpha1 = obr_.lookupObject<volScalarField>("alpha1");

const fvMesh& mesh = U.mesh();
     
        const surfaceVectorField::GeometricBoundaryField& Sfb =
            mesh.Sf().boundaryField();
            
            
        tmp<volSymmTensorField> tdevRhoReff = devRhoReff();
     const volSymmTensorField::GeometricBoundaryField& devRhoReffb
            = tdevRhoReff().boundaryField();
                        
        const volScalarField limitedAlpha1
        (
           min(max(alpha1, scalar(0)), scalar(1))
        );

forAllConstIter(labelHashSet, patchSet_, iter)
        {   
            label patchi = iter.key();

            
vectorField pf1(Sfb[patchi]*(p.boundaryField([patchi]*limitedAlpha1.boundaryField()[patchi]));
vectorField pf2(Sfb[patchi]*(p.boundaryField()[patchi]*(scalar(1)-limitedAlpha1.boundaryField()[patchi])));


vectorField vf1((Sfb[patchi]*limitedAlpha1.boundaryField()[patchi]) & devRhoReffb[patchi]);
vectorField vf2((Sfb[patchi]*(scalar(1)-limitedAlpha1.boundaryField()[patchi])) & devRhoReffb[patchi]);



fmPressurePhase1 += sum(pf1);
fmPressurePhase2 += sum(pf2);


fmViscousPhase1 += sum(vf1);
fmViscousPhase2 += sum(vf2);

and some other changes to print them out. Don't know if it makes sense, i need to test it. I think a standard test case could be poiseuille flow between two parallel plates with 2 fluid and a meniscus between them. When there is only phase 1 (at the beginning) or 2(at the end of sim) i should find analytical solution for viscous dissipation.
What do you think?

best

andrea