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

Force calculation in multiphase simulations

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

Like Tree8Likes

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   March 1, 2012, 05:11
Default Force calculation in multiphase simulations
  #1
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 319
Rep Power: 16
Andrea_85 is on a distinguished road
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
keepfit likes this.
Andrea_85 is offline   Reply With Quote

Old   March 2, 2012, 04:20
Default
  #2
Senior Member
 
niaz's Avatar
 
A_R
Join Date: Jun 2009
Posts: 122
Rep Power: 16
niaz is on a distinguished road
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.
niaz is offline   Reply With Quote

Old   March 2, 2012, 05:24
Default
  #3
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 319
Rep Power: 16
Andrea_85 is on a distinguished road
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
Andrea_85 is offline   Reply With Quote

Old   March 2, 2012, 07:32
Default
  #4
Senior Member
 
calim_cfd's Avatar
 
mauricio
Join Date: Jun 2011
Posts: 172
Rep Power: 17
calim_cfd is on a distinguished road
Quote:
Originally Posted by Andrea_85 View Post
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
calim_cfd is offline   Reply With Quote

Old   March 3, 2012, 02:53
Default
  #5
Senior Member
 
niaz's Avatar
 
A_R
Join Date: Jun 2009
Posts: 122
Rep Power: 16
niaz is on a distinguished road
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.
niaz is offline   Reply With Quote

Old   March 5, 2012, 12:25
Default
  #6
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 319
Rep Power: 16
Andrea_85 is on a distinguished road
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
waku2005, dupeng and Tolga KURUMUS like this.
Andrea_85 is offline   Reply With Quote

Old   March 7, 2012, 20:17
Default
  #7
Senior Member
 
Dave
Join Date: Jul 2010
Posts: 100
Rep Power: 15
daveatstyacht is on a distinguished road
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
daveatstyacht is offline   Reply With Quote

Old   March 8, 2012, 05:03
Default
  #8
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 319
Rep Power: 16
Andrea_85 is on a distinguished road
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
Andrea_85 is offline   Reply With Quote

Old   March 8, 2012, 11:22
Default
  #9
Senior Member
 
Dave
Join Date: Jul 2010
Posts: 100
Rep Power: 15
daveatstyacht is on a distinguished road
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
daveatstyacht is offline   Reply With Quote

Old   March 16, 2012, 10:44
Default
  #10
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 319
Rep Power: 16
Andrea_85 is on a distinguished road
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
Attached Images
File Type: jpg forces.jpg (44.8 KB, 377 views)
Andrea_85 is offline   Reply With Quote

Old   March 16, 2012, 14:01
Default
  #11
Senior Member
 
calim_cfd's Avatar
 
mauricio
Join Date: Jun 2011
Posts: 172
Rep Power: 17
calim_cfd is on a distinguished road
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
__________________
Best Regards
/calim

"Elune will grant us the strength"
calim_cfd is offline   Reply With Quote

Old   March 19, 2012, 12:27
Default
  #12
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 319
Rep Power: 16
Andrea_85 is on a distinguished road
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
Andrea_85 is offline   Reply With Quote

Old   March 19, 2012, 13:25
Default
  #13
Senior Member
 
calim_cfd's Avatar
 
mauricio
Join Date: Jun 2011
Posts: 172
Rep Power: 17
calim_cfd is on a distinguished road
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
__________________
Best Regards
/calim

"Elune will grant us the strength"
calim_cfd is offline   Reply With Quote

Old   March 21, 2012, 10:44
Default
  #14
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 319
Rep Power: 16
Andrea_85 is on a distinguished road
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
Attached Files
File Type: zip fallingSlug.zip (15.3 KB, 96 views)
Andrea_85 is offline   Reply With Quote

Old   March 21, 2012, 10:46
Default
  #15
Senior Member
 
calim_cfd's Avatar
 
mauricio
Join Date: Jun 2011
Posts: 172
Rep Power: 17
calim_cfd is on a distinguished road
hi andrea!

thx 4 sharing!
i'll give it a try and let you know if sth comes up
__________________
Best Regards
/calim

"Elune will grant us the strength"
calim_cfd is offline   Reply With Quote

Old   March 21, 2012, 12:35
Default
  #16
Senior Member
 
calim_cfd's Avatar
 
mauricio
Join Date: Jun 2011
Posts: 172
Rep Power: 17
calim_cfd is on a distinguished road
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

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
__________________
Best Regards
/calim

"Elune will grant us the strength"

Last edited by calim_cfd; March 22, 2012 at 06:28.
calim_cfd is offline   Reply With Quote

Old   March 22, 2012, 12:12
Default
  #17
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 319
Rep Power: 16
Andrea_85 is on a distinguished road
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
Attached Images
File Type: jpg thinFilm.jpg (44.9 KB, 139 views)
Andrea_85 is offline   Reply With Quote

Old   March 22, 2012, 12:48
Default
  #18
Senior Member
 
calim_cfd's Avatar
 
mauricio
Join Date: Jun 2011
Posts: 172
Rep Power: 17
calim_cfd is on a distinguished road
Quote:
Originally Posted by Andrea_85 View Post
....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!
__________________
Best Regards
/calim

"Elune will grant us the strength"
calim_cfd is offline   Reply With Quote

Old   March 28, 2012, 09:09
Default
  #19
Member
 
Björn Windén
Join Date: Feb 2012
Location: National Maritime Research Institute, Tokyo, Japan
Posts: 37
Rep Power: 14
winden is on a distinguished road
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
Attached Images
File Type: jpg anim1.0007.jpg (54.5 KB, 186 views)
winden is offline   Reply With Quote

Old   June 4, 2012, 09:33
Default
  #20
Senior Member
 
niaz's Avatar
 
A_R
Join Date: Jun 2009
Posts: 122
Rep Power: 16
niaz is on a distinguished road
Dear Andrea
can you send your library for me?

Quote:
Originally Posted by Andrea_85 View Post
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
niaz 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
Calculation of Drag Force per unit Length Mohammad Faridul Alam CFX 4 January 11, 2013 08:19
Time-averaged force in CFX rjmcsherry CFX 2 October 21, 2010 11:34
Calculation of Saffman Lift Force JPBodner Main CFD Forum 3 August 4, 2010 12:11
MRF and Heat transfer calculation Susan YU FLUENT 0 June 2, 2010 09:46
Viscous Force Calculation cwang5 OpenFOAM Programming & Development 1 May 4, 2010 05:59


All times are GMT -4. The time now is 00:38.