# Force calculation in multiphase simulations

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

 March 1, 2012, 05:11 Force calculation in multiphase simulations #1 Senior Member   Andrea Ferrari Join Date: Dec 2010 Posts: 319 Rep Power: 12 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.

 March 2, 2012, 04:20 #2 Senior Member     A_R Join Date: Jun 2009 Posts: 122 Rep Power: 13 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.

 March 2, 2012, 05:24 #3 Senior Member   Andrea Ferrari Join Date: Dec 2010 Posts: 319 Rep Power: 12 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

March 2, 2012, 07:32
#4
Senior Member

mauricio
Join Date: Jun 2011
Posts: 151
Rep Power: 14
Quote:
 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

 March 3, 2012, 02:53 #5 Senior Member     A_R Join Date: Jun 2009 Posts: 122 Rep Power: 13 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.

 March 5, 2012, 12:25 #6 Senior Member   Andrea Ferrari Join Date: Dec 2010 Posts: 319 Rep Power: 12 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 waku2005, dupeng and Tolga KURUMUS like this.

 March 7, 2012, 20:17 #7 Member   Dave Join Date: Jul 2010 Posts: 98 Rep Power: 12 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

 March 8, 2012, 05:03 #8 Senior Member   Andrea Ferrari Join Date: Dec 2010 Posts: 319 Rep Power: 12 Hi Dave, you are right but i'm working with multiphase flow and my cases are always laminar (0.01

 March 8, 2012, 11:22 #9 Member   Dave Join Date: Jul 2010 Posts: 98 Rep Power: 12 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

March 16, 2012, 10:44
#10
Senior Member

Andrea Ferrari
Join Date: Dec 2010
Posts: 319
Rep Power: 12
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.

best regards

andrea
Attached Images
 forces.jpg (44.8 KB, 333 views)

 March 16, 2012, 14:01 #11 Senior Member     mauricio Join Date: Jun 2011 Posts: 151 Rep Power: 14 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"

 March 19, 2012, 12:27 #12 Senior Member   Andrea Ferrari Join Date: Dec 2010 Posts: 319 Rep Power: 12 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

 March 19, 2012, 13:25 #13 Senior Member     mauricio Join Date: Jun 2011 Posts: 151 Rep Power: 14 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"

March 21, 2012, 10:44
#14
Senior Member

Andrea Ferrari
Join Date: Dec 2010
Posts: 319
Rep Power: 12
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
 fallingSlug.zip (15.3 KB, 91 views)

 March 21, 2012, 10:46 #15 Senior Member     mauricio Join Date: Jun 2011 Posts: 151 Rep Power: 14 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"

 March 21, 2012, 12:35 #16 Senior Member     mauricio Join Date: Jun 2011 Posts: 151 Rep Power: 14 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.

March 22, 2012, 12:12
#17
Senior Member

Andrea Ferrari
Join Date: Dec 2010
Posts: 319
Rep Power: 12
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
 thinFilm.jpg (44.9 KB, 126 views)

March 22, 2012, 12:48
#18
Senior Member

mauricio
Join Date: Jun 2011
Posts: 151
Rep Power: 14
Quote:
 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!
__________________
Best Regards
/calim

"Elune will grant us the strength"

March 28, 2012, 09:09
#19
Member

Björn Windén
Join Date: Feb 2012
Location: National Maritime Research Institute, Tokyo, Japan
Posts: 37
Rep Power: 10
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
 anim1.0007.jpg (54.5 KB, 173 views)

June 4, 2012, 09:33
#20
Senior Member

A_R
Join Date: Jun 2009
Posts: 122
Rep Power: 13
Dear Andrea
can you send your library for me?

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