
[Sponsors] 
March 1, 2012, 05:11 
Force calculation in multiphase simulations

#1 
Senior Member
Andrea Ferrari
Join Date: Dec 2010
Posts: 296
Rep Power: 8 
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/OpenFOAM2.1.0/tutorials/multiphase/interFoam/laminar/damBreak/constant/transportProperties" file: /home/aferrari/OpenFOAM/OpenFOAM2.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 

March 2, 2012, 04:20 

#2 
Senior Member
A_R
Join Date: Jun 2009
Posts: 118
Rep Power: 9 
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: 296
Rep Power: 8 
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: 141
Rep Power: 9 
Quote:
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: 118
Rep Power: 9 
Andrea
if you look at base file for force and forcecoeff you can find that both are for singlephase 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: 296
Rep Power: 8 
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))); 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); What do you think? best andrea 

March 7, 2012, 20:17 

#7 
Member
Dave
Join Date: Jul 2010
Posts: 97
Rep Power: 8 
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 underpredict forces for flows with large mixture regions. Dave 

March 8, 2012, 05:03 

#8 
Senior Member
Andrea Ferrari
Join Date: Dec 2010
Posts: 296
Rep Power: 8 
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 

March 8, 2012, 11:22 

#9 
Member
Dave
Join Date: Jul 2010
Posts: 97
Rep Power: 8 
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: 296
Rep Power: 8 
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 

March 16, 2012, 14:01 

#11 
Senior Member
mauricio
Join Date: Jun 2011
Posts: 141
Rep Power: 9 
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: 296
Rep Power: 8 
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: 141
Rep Power: 9 
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=qGRnT4fBK8iN4gSj_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: 296
Rep Power: 8 
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 

March 21, 2012, 10:46 

#15 
Senior Member
mauricio
Join Date: Jun 2011
Posts: 141
Rep Power: 9 
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: 141
Rep Power: 9 
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; } 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: 296
Rep Power: 8 
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 

March 22, 2012, 12:48 

#18  
Senior Member
mauricio
Join Date: Jun 2011
Posts: 141
Rep Power: 9 
Quote:
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, 08:09 

#19 
Member
Björn Windén
Join Date: Feb 2012
Location: University of Southampton UK
Posts: 37
Rep Power: 6 
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 phasedictionaries 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 forcelibrary 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 

June 4, 2012, 08:33 

#20  
Senior Member
A_R
Join Date: Jun 2009
Posts: 118
Rep Power: 9 
Dear Andrea
can you send your library for me? Quote:


Thread Tools  
Display Modes  


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 
Timeaveraged force in CFX  rjmcsherry  CFX  2  October 21, 2010 10:34 
Calculation of Saffman Lift Force  JPBodner  Main CFD Forum  3  August 4, 2010 11:11 
MRF and Heat transfer calculation  Susan YU  FLUENT  0  June 2, 2010 08:46 
Viscous Force Calculation  cwang5  OpenFOAM Programming & Development  1  May 4, 2010 04:59 