
[Sponsors] 
September 28, 2009, 20:04 
Calculate Lift and Drag Coefficients CL and CD

#1 
Member
Sven Winkler
Join Date: May 2009
Posts: 70
Rep Power: 9 
Hi,
I am simulating a turbulent flow around a square. For This flow, I want to calculate the lift and the drag coefficients CD and CL. By scanning through posts in this forum, I found that this is possible by adding the following lines to the controlDict file: Code:
functions ( forces { type forces; functionObjectLibs ("libforces.so"); //Lib to load patches (cylinder); // change to your patch name rhoInf 1.204; //Reference density for fluid CofR (0 0 0); //Origin for moment calculations outputControl timeStep; outputInterval 100; } forceCoeffs { type forceCoeffs; functionObjectLibs ("libforces.so"); patches (cylinder); //change to your patch name rhoInf 1.204; CofR (0 0 0); liftDir (0 1 0); dragDir (1 0 0); pitchAxis (0 0 0); magUInf 0.1; lRef 1; Aref 1; outputControl timeStep; outputInterval 100; } );


September 28, 2009, 22:17 

#2 
Member
Simon Lapointe
Join Date: May 2009
Location: Québec, Qc, Canada
Posts: 33
Rep Power: 9 
Hi,
1. rhoInf is the freestream density. 2. magUinf is the magnitude of the freestream velocity. 3. The source code for the forces and forceCoeffs are located, respectively, in OpenFOAM1.6/src/postProcessing/forces/ and OpenFOAM1.6/src/postProcessing/forceCoeffs/ 4. The equations of CL and CD are defined in OpenFOAM1.6/src/postProcessing/forceCoeffs/forceCoeffs.C Hope this helps 

September 28, 2009, 22:54 

#3 
Member
Sven Winkler
Join Date: May 2009
Posts: 70
Rep Power: 9 
Thanks for your answer Simon. I still have a problem with magUInf, the magnitude of the free stream velocity. What exactly is the magnitude of my free stream velocity? My free stream velocity is 0.322, so the magnitude is 0.1? I tried it with that value (0.1), but it does not yield the correct results! Thanks again for your help!


September 29, 2009, 02:53 

#4 
New Member
Matthew Philpott
Join Date: Aug 2009
Location: Belgium
Posts: 24
Rep Power: 9 
If your free stream velocity is 0.322, then I'm guessing the magnitude is 0.322. If that is the only velocity component ie, value (0.322 0 0). Magnitude is calculated using typical vector mathematics v²=vx²+vy²+vz². The square of the velocity magnitude is the sum of the squares of your velocity components.
__________________
CAELinux 2009 + OF1.5 Ubuntu 9.04 x64 (jaunty jackalope) + OF1.6 

September 29, 2009, 21:01 

#5 
Member
Sven Winkler
Join Date: May 2009
Posts: 70
Rep Power: 9 
Ok..now I got it! Thank you all very much for your help. Greetings from California
Sven 

October 2, 2009, 21:19 

#6 
Member
Sven Winkler
Join Date: May 2009
Posts: 70
Rep Power: 9 
Well, I looked up the definition for CD and CL in the source Code. First of all, the right location for the files is:
/home/user/OpenFOAM/OpenFOAM1.6/src/postProcessing/functionObjects/forces There are two folders (forces and forceCoeffs) where you can find the appropriate .C and .H files. I first had a look at the forceCoeffs.C file, here is part of that: Code:
void Foam::forceCoeffs::write() { if (active_) { // Create the forces file if not already created makeFile(); forcesMoments fm = forces::calcForcesMoment(); scalar pDyn = 0.5*rhoRef_*magUInf_*magUInf_; vector totForce = fm.first().first() + fm.first().second(); vector totMoment = fm.second().first() + fm.second().second(); scalar liftForce = totForce & liftDir_; scalar dragForce = totForce & dragDir_; scalar pitchMoment = totMoment & pitchAxis_; scalar Cl = liftForce/(Aref_*pDyn); scalar Cd = dragForce/(Aref_*pDyn); scalar Cm = pitchMoment/(Aref_*lRef_*pDyn); if (Pstream::master()) { forcesFilePtr_() << obr_.time().value() << tab << Cd << tab << Cl << tab << Cm << endl; if (log_) { Info<< "forceCoeffs output:" << nl << " Cd = " << Cd << nl << " Cl = " << Cl << nl << " Cm = " << Cm << nl << endl; } } } } Unfortunately I do not understand, how OpenFOAM calculates the total force mentioned above. I tried to look up the definition of first() and second(), which are defined in Tuple2.H (OpenFoam/InInclude/Tuple2.H). Here is the definition of these functions: Code:
// Return first inline const Type1& first() const { return f_; } // Return first inline Type1& first() { return f_; } // Return second inline const Type2& second() const { return s_; } // Return second inline Type2& second() { return s_; } 

October 2, 2009, 21:21 

#7 
Member
Sven Winkler
Join Date: May 2009
Posts: 70
Rep Power: 9 
In the forces.C file, I also found this:
Code:
Foam::forces::forcesMoments Foam::forces::calcForcesMoment() const { forcesMoments fm ( pressureViscous(vector::zero, vector::zero), pressureViscous(vector::zero, vector::zero) ); if (directForceDensity_) { const volVectorField& fD = obr_.lookupObject<volVectorField>(fDName_); const fvMesh& mesh = fD.mesh(); const surfaceVectorField::GeometricBoundaryField& Sfb = mesh.Sf().boundaryField(); forAllConstIter(labelHashSet, patchSet_, iter) { label patchi = iter.key(); vectorField Md = mesh.C().boundaryField()[patchi]  CofR_; scalarField sA = mag(Sfb[patchi]); // Normal force = surfaceUnitNormal * (surfaceNormal & forceDensity) vectorField fN = Sfb[patchi]/sA *( Sfb[patchi] & fD.boundaryField()[patchi] ); fm.first().first() += sum(fN); fm.second().first() += sum(Md ^ fN); // Tangential force (total force minus normal fN) vectorField fT = sA*fD.boundaryField()[patchi]  fN; fm.first().second() += sum(fT); fm.second().second() += sum(Md ^ fT); } } else { const volVectorField& U = obr_.lookupObject<volVectorField>(UName_); const volScalarField& p = obr_.lookupObject<volScalarField>(pName_); const fvMesh& mesh = U.mesh(); const surfaceVectorField::GeometricBoundaryField& Sfb = mesh.Sf().boundaryField(); tmp<volSymmTensorField> tdevRhoReff = devRhoReff(); const volSymmTensorField::GeometricBoundaryField& devRhoReffb = tdevRhoReff().boundaryField(); forAllConstIter(labelHashSet, patchSet_, iter) { label patchi = iter.key(); vectorField Md = mesh.C().boundaryField()[patchi]  CofR_; vectorField pf = Sfb[patchi]*p.boundaryField()[patchi]; fm.first().first() += rho(p)*sum(pf); fm.second().first() += rho(p)*sum(Md ^ pf); vectorField vf = Sfb[patchi] & devRhoReffb[patchi]; fm.first().second() += sum(vf); fm.second().second() += sum(Md ^ vf); } } reduce(fm, sumOp()); return fm; } Thanks a lot for all your help! 

August 10, 2011, 05:21 

#8  
Member

Quote:
_____ Rasoul 

June 13, 2012, 06:20 

#9 
New Member
Malhar Malushte
Join Date: May 2012
Posts: 16
Rep Power: 6 
hello every1,
i am also new to openfoam. i have done similar sim ulation for flow acrodss a cylinder. i get drag n lift forces correctly, i.e accor ding to vortex shedding i am getting variation in lift forces. but the coeff of drag n lift Cd n Cl, i am getting them constant throughout. i am unable ti digest this contrasting behavour. please guide me thr this. thanks n regards malhar. 

December 17, 2012, 05:53 

#10  
Member
Malik
Join Date: Dec 2012
Location: Austin, USA
Posts: 52
Rep Power: 6 
Quote:
by saying this, do you mean that the pressure value used is the one on the boundary of the cells which touch the surface on which drag and lift force are applied ? I am a new Foamer but i think that we forget the viscous force if OpenFOAM only does this. In the code, there is a tangential force calculated, but i cannot say if it concerns the viscous force or not. Does anyone have an idea ? 

August 1, 2014, 01:29 

#11 
Senior Member

Understanding of the terms...
I am actually looking out for the term mentioned below in the previous posts on this thread: const volVectorField& fD = obr_.lookupObject<volVectorField>(fDName_); Where is the term located in the case directory? Correct me if I am wrong Thanks in Advance! 

August 1, 2014, 01:38 

#12 
Senior Member

Oh, it's clear the term is for compressible flows.
Please, excuse for previous post. I did mistake it for incompressible flows. Last edited by Tushar@cfd; August 2, 2014 at 01:24. 

Thread Tools  
Display Modes  

