CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)
-   -   Force calculation in multiphase simulations (https://www.cfd-online.com/Forums/openfoam-solving/98003-force-calculation-multiphase-simulations.html)

 Andrea_85 March 1, 2012 05:11

Force calculation in multiphase simulations

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 March 2, 2012 04:20

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 March 2, 2012 05:24

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 March 2, 2012 07:32

Quote:
 Originally Posted by Andrea_85 (Post 347247) 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 March 3, 2012 02:53

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 March 5, 2012 12:25

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

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

Code:

```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()

Code:

```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 March 7, 2012 20:17

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 March 8, 2012 05:03

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 March 8, 2012 11:22

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 March 16, 2012 10:44

1 Attachment(s)
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..:confused:). 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.

best regards

andrea

 calim_cfd March 16, 2012 14:01

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

 Andrea_85 March 19, 2012 12:27

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.:confused:

best

andrea

 calim_cfd March 19, 2012 13:25

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 March 21, 2012 10:44

1 Attachment(s)
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 March 21, 2012 10:46

hi andrea!

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

 calim_cfd March 21, 2012 12:35

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

``` 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
http://s13.postimage.org/5r93roi03/error.jpg
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 March 22, 2012 12:12

1 Attachment(s)
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 March 22, 2012 12:48

Quote:
 Originally Posted by Andrea_85 (Post 350944) ....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 (?) :D

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 March 28, 2012 08:09

1 Attachment(s)
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 June 4, 2012 08:33

Dear Andrea
can you send your library for me?

Quote:
 Originally Posted by Andrea_85 (Post 347744) Hi Niaz and Calim, thanks for your suggestions. This is what i've done: in forces.C where devRhoReff() is calculated i added: Code: ```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() Code: ```const volVectorField& U = obr_.lookupObject(UName_); const volScalarField& p = obr_.lookupObject(pName_); const volScalarField& alpha1 = obr_.lookupObject("alpha1"); const fvMesh& mesh = U.mesh();             const surfaceVectorField::GeometricBoundaryField& Sfb =             mesh.Sf().boundaryField();                                 tmp 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

All times are GMT -4. The time now is 11:57.