# Drag coefficient implementation

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

 November 30, 2020, 02:44 Drag coefficient implementation #1 New Member   Sumit Peh Join Date: Oct 2018 Location: Beijing Posts: 20 Rep Power: 7 Hi Foamers I would like to apply equations of other Drag coefficient in 4 conditions of Renold's number. [Re<1 // 11000] SchillerNaumann has 2 conditions of Re which are Re<1000 and Re>1000 Here is the SchillerNaumann code on ~dragModels I supposed neg(Re - 1000) is this Re<1000 condition and pos0(Re - 1000) is the Re>1000. Code: { volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3))); volScalarField Cds ( neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re) + pos0(Re - 1000)*0.44 ); return 0.75*Cds*phase2_.rho()*Ur/phase1_.d(); } Is there a way to write code for more than 2 conditions ? As my imagination, I guess maybe my case would be Code: { volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3))); volScalarField Cds ( neg(Re - 1)*(Equation1) + neg1(Re - 100)*(Equation2) + neg100(Re - 1000)*(Equation3) + pos0(Re - 1000)*(Equation4) ); return 0.75*Cds*phase2_.rho()*Ur/phase1_.d(); } Please confirm if i'm right and please correct if i'm wrong. Thanks in advanced James

 November 30, 2020, 04:04 #2 Senior Member   Join Date: Dec 2019 Location: Cologne, Germany Posts: 358 Rep Power: 8 if the compiler does not complain about your code, you still need to verify it with easy to check simulation, maybe a quadratic pipe with one cell depth and height and 30 cells long. set different inlet velocities for both phases and use the calculator in paraview to calculate your drag coeff., also you can output the drag coeff. with openfoam, so you can cross check.

December 7, 2020, 22:02
#3
New Member

Sumit Peh
Join Date: Oct 2018
Location: Beijing
Posts: 20
Rep Power: 7
Quote:
 Originally Posted by geth03 if the compiler does not complain about your code, you still need to verify it with easy to check simulation, maybe a quadratic pipe with one cell depth and height and 30 cells long. set different inlet velocities for both phases and use the calculator in paraview to calculate your drag coeff., also you can output the drag coeff. with openfoam, so you can cross check.
So i tried this instead. Everything go well when compile but not sure its correct or not.

+ (1 - Re -100)
+ (100 - Re - 1000)

 December 9, 2020, 03:52 #4 Senior Member   Join Date: Dec 2019 Location: Cologne, Germany Posts: 358 Rep Power: 8 ok, let us first clarify what neg() and pos() mean: neg(scalar s) checks if the value of s i negative, if yes, it will return the value 1. so if neg(Re-1000) is negativ, than you will get a 1 back. the same goes for pos(s) or pos0(s). pos(s): is s larger than zero, pos0(s): is s larger than zero or 0. if yes, return 1. which means, if you have neg(s) and pos0(s) in your code, it will only return for one case 1, otherwise 0. if you have many of that, as you did, it can and will be the case, that you will return more 1 than expected. let me give an example: if Re = 10 and you have + neg(1 - Re -100) + neg(100 - Re - 1000) like you did, you will have +1 +1 as return values, because in both cases, the if-check will give "true". so, in that case, you will overpredict the drag coefficient. as a workaround, i suggest you loop over all your cells and depending on the Re, you can check in which value area Re is for each cell and save Cds accordingly. afterwards you can insert that one in the place where this short formula would normally be. i hope things are clear now.

December 10, 2020, 03:50
#5
New Member

Sumit Peh
Join Date: Oct 2018
Location: Beijing
Posts: 20
Rep Power: 7
Quote:
 Originally Posted by geth03 ok, let us first clarify what neg() and pos() mean: neg(scalar s) checks if the value of s i negative, if yes, it will return the value 1. so if neg(Re-1000) is negativ, than you will get a 1 back. the same goes for pos(s) or pos0(s). pos(s): is s larger than zero, pos0(s): is s larger than zero or 0. if yes, return 1. which means, if you have neg(s) and pos0(s) in your code, it will only return for one case 1, otherwise 0. if you have many of that, as you did, it can and will be the case, that you will return more 1 than expected. let me give an example: if Re = 10 and you have + neg(1 - Re -100) + neg(100 - Re - 1000) like you did, you will have +1 +1 as return values, because in both cases, the if-check will give "true". so, in that case, you will overpredict the drag coefficient. as a workaround, i suggest you loop over all your cells and depending on the Re, you can check in which value area Re is for each cell and save Cds accordingly. afterwards you can insert that one in the place where this short formula would normally be. i hope things are clear now.
Hi geth03,

I also found a Lain.C in /drags, it use the code like this so I changed it.

Code:
```#include "Lain.H"
#include "phasePair.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
namespace dragModels
{
defineTypeNameAndDebug(Lain, 0);
}
}

// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::dragModels::Lain::Lain
(
const dictionary& dict,
const phasePair& pair,
const bool registerObject
)
:
dragModel(dict, pair, registerObject)
{}

// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::dragModels::Lain::~Lain()
{}

// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

Foam::tmp<Foam::volScalarField> Foam::dragModels::Lain::CdRe() const
{
volScalarField Re(pair_.Re());

return
neg(Re - 1.5)*16.0
+ pos0(Re - 1.5)*neg(Re - 80)*14.9*pow(Re, 0.22)
+ pos0(Re - 80)*neg(Re - 1500)*48*(1 - 2.21/sqrt(max(Re, SMALL)))
+ pos0(Re - 1500)*2.61*Re;
}```
As your example, why Re = 10 get +1 result ? It isn't in the range of neg(100 - Re - 1000).

According to my understanding, we can write code by this way :
pos0(Re - 1)*neg(Re - 100)
for a condition of 1 < Re < 100

 December 10, 2020, 05:28 #6 Senior Member   Join Date: Dec 2019 Location: Cologne, Germany Posts: 358 Rep Power: 8 if Re = 10, then neg(100 - Re - 1000) = neg(-910). neg(-910) will check if -910 is a negative value, and return +1 if true, which it is, in that case. as in the Lain.C - case, you can use pos0(s) neg(s1) in combination so if your Re falls between two transition values, both return +1 and +1, otherwise +1 and 0, which makes +1 x 0 = 0, and so other Re-ranges are neglected. i didn't think about combining those pos() and neg(), that is way smarter then looping over all cells. Jamessmp23 likes this.

 Tags cfd, code, drag, drag coeffcients, openfoam