Permalink
Browse files

twoPhaseEulerFoam: Corrected the drag models

such that the drag-coefficient is correct on processor patches

Implements patches provided and closes bug-report
http://www.openfoam.com/mantisbt/view.php?id=205
  • Loading branch information...
1 parent ab87047 commit 1e7d927e96835885f4824e750bd182f76424f458 Henry committed May 25, 2011
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
+ \\ / A nd | Copyright (C) 1991-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -73,34 +73,19 @@ Foam::tmp<Foam::volScalarField> Foam::GidaspowErgunWenYu::K
volScalarField bp = pow(beta, -2.65);
volScalarField Re = max(Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3));
- volScalarField Cds = 24.0*(1.0 + 0.15*pow(Re, 0.687))/Re;
+ volScalarField Cds =
+ neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+ + pos(Re - 1000)*0.44;
- forAll(Re, celli)
- {
- if(Re[celli] > 1000.0)
- {
- Cds[celli] = 0.44;
- }
- }
-
// Wen and Yu (1966)
- tmp<volScalarField> tKWenYu = 0.75*Cds*phaseb_.rho()*Ur*bp/phasea_.d();
- volScalarField& KWenYu = tKWenYu();
-
- // Ergun
- forAll (beta, cellj)
- {
- if (beta[cellj] <= 0.8)
- {
- KWenYu[cellj] =
- 150.0*alpha_[cellj]*phaseb_.nu().value()*phaseb_.rho().value()
- /sqr(beta[cellj]*phasea_.d().value())
- + 1.75*phaseb_.rho().value()*Ur[cellj]
- /(beta[cellj]*phasea_.d().value());
- }
- }
-
- return tKWenYu;
+ return
+ pos(beta - 0.8)
+ *(0.75*alpha_*beta*Cds*phaseb_.rho()*Ur*bp/phasea_.d())
+ + neg(beta - 0.8)
+ *(
+ 150.0*sqr(alpha_)*phaseb_.nu()*phaseb_.rho()/(beta*sqr(phasea_.d()))
+ + 1.75*alpha_*phaseb_.rho()*Ur/phasea_.d()
+ );
}
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
+ \\ / A nd | Copyright (C) 1991-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -72,15 +72,10 @@ Foam::tmp<Foam::volScalarField> Foam::GidaspowSchillerNaumann::K
volScalarField bp = pow(beta, -2.65);
volScalarField Re = max(beta*Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3));
- volScalarField Cds = 24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re;
-
- forAll(Re, celli)
- {
- if(Re[celli] > 1000.0)
- {
- Cds[celli] = 0.44;
- }
- }
+
+ volScalarField Cds =
+ neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+ + pos(Re - 1000)*0.44;
return 0.75*Cds*phaseb_.rho()*Ur*bp/phasea_.d();
}
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
+ \\ / A nd | Copyright (C) 1991-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -69,15 +69,10 @@ Foam::tmp<Foam::volScalarField> Foam::SchillerNaumann::K
) const
{
volScalarField Re = max(Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3));
- volScalarField Cds = 24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re;
-
- forAll(Re, celli)
- {
- if(Re[celli] > 1000.0)
- {
- Cds[celli] = 0.44;
- }
- }
+
+ volScalarField Cds =
+ neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+ + pos(Re - 1000)*0.44;
return 0.75*Cds*phaseb_.rho()*Ur/phasea_.d();
}
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
+ \\ / A nd | Copyright (C) 1991-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -70,15 +70,9 @@ Foam::tmp<Foam::volScalarField> Foam::SyamlalOBrien::K
{
volScalarField beta = max(scalar(1) - alpha_, scalar(1.0e-6));
volScalarField A = pow(beta, 4.14);
- volScalarField B = 0.8*pow(beta, 1.28);
-
- forAll (beta, celli)
- {
- if (beta[celli] > 0.85)
- {
- B[celli] = pow(beta[celli], 2.65);
- }
- }
+ volScalarField B =
+ neg(beta - 0.85)*(0.8*pow(beta, 1.28))
+ + pos(beta - 0.85)*(pow(beta, 2.65));
volScalarField Re = max(Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3));
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
+ \\ / A nd | Copyright (C) 1991-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -72,15 +72,9 @@ Foam::tmp<Foam::volScalarField> Foam::WenYu::K
volScalarField bp = pow(beta, -2.65);
volScalarField Re = max(Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3));
- volScalarField Cds = 24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re;
-
- forAll(Re, celli)
- {
- if(Re[celli] > 1000.0)
- {
- Cds[celli] = 0.44;
- }
- }
+ volScalarField Cds =
+ neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+ + pos(Re - 1000)*0.44;
return 0.75*Cds*phaseb_.rho()*Ur*bp/phasea_.d();
}

0 comments on commit 1e7d927

Please sign in to comment.