View Issue Details

Jump to NotesJump to History
IDProjectCategoryView StatusDate SubmittedLast Update
0000205OpenFOAMBugpublic2011-05-25 00:352011-05-26 08:27
Reporteralbertop Assigned Tohenry  
PrioritynormalSeveritymajorReproducibilityalways
Status resolvedResolutionfixed 
PlatformLinuxOSOpenSuseOS Version11.3
Summary0000205: Inconsistent behavior at processor boundaries in twoPhaseEulerFoam
DescriptionThe twoPhaseEulerFoam solver shows inconsistencies immediately after the processor boundaries used in parallel cases. This translates in altered values for the fields (alpha, U, ...) after such values.
Steps To ReproduceRun twoPhaseEulerFoam on the testcase bed2 on 4 processors, with simple decomposition along the y direction. Visualize the field of alpha in the portion of the system with 0.2 < y < 0.4, and rescale alpha between 0.54 and 0.55 at t = 0.01.
The case is attached for your convenience, with the decomposeParDict.
Additional InformationThe tar.gz file contains two pictures, showing the results in serial and in parallel.
TagsNo tags attached.

Activities

albertop

2011-05-25 00:35

reporter  

bed2.tar.gz (31,624 bytes)

albertop

2011-05-25 09:57

reporter  

GidaspowErgunWenYu.patch (1,586 bytes)    
*** GidaspowErgunWenYu.C	2011-05-25 03:55:26.696938481 -0500
--- GidaspowErgunWenYuPatched.C	2011-05-25 03:55:18.764643670 -0500
***************
*** 73,106 ****
      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;
! 
!     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;
  }
  
  
--- 73,86 ----
      volScalarField bp = pow(beta, -2.65);
      volScalarField Re = max(Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3));
  
!     volScalarField Cds = 
! 	neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
!       + pos(Re - 1000)*0.44;
!  
      // Wen and Yu (1966)
!     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());
  }
GidaspowErgunWenYu.patch (1,586 bytes)   

albertop

2011-05-25 09:58

reporter   ~0000392

It seems the problem was caused by the implementation of the "GidaspowErgunWenYu" drag model, which did not update value properly at boundaries.

I attach a patch. I will look into the other drag models and let you know.

Thank you.

albertop

2011-05-25 10:43

reporter   ~0000393

Similar patches for Cds calculation to remove the loop in

- interfacialModel/GidaspowSchillerNaumann
- interfacialModel/SchillerNaumann
- interfacialModel/WenYu

and for B in interfacialModel/SyamlalOBrien.

I directly attach the corrected .C files.

albertop

2011-05-25 10:44

reporter  

GidaspowSchillerNaumannPatched.C (2,601 bytes)    
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "GidaspowSchillerNaumann.H"
#include "addToRunTimeSelectionTable.H"

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

namespace Foam
{
    defineTypeNameAndDebug(GidaspowSchillerNaumann, 0);

    addToRunTimeSelectionTable
    (
        dragModel,
        GidaspowSchillerNaumann,
        dictionary
    );
}


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

Foam::GidaspowSchillerNaumann::GidaspowSchillerNaumann
(
    const dictionary& interfaceDict,
    const volScalarField& alpha,
    const phaseModel& phasea,
    const phaseModel& phaseb
)
:
    dragModel(interfaceDict, alpha, phasea, phaseb)
{}


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

Foam::GidaspowSchillerNaumann::~GidaspowSchillerNaumann()
{}


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

Foam::tmp<Foam::volScalarField> Foam::GidaspowSchillerNaumann::K
(
    const volScalarField& Ur
) const
{
    volScalarField beta = max(scalar(1) - alpha_, scalar(1e-6));
    volScalarField bp = pow(beta, -2.65);

    volScalarField Re = max(beta*Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3));
    
    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();
}


// ************************************************************************* //
GidaspowSchillerNaumannPatched.C (2,601 bytes)   

albertop

2011-05-25 10:44

reporter  

SchillerNaumannPatched.C (2,421 bytes)    
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "SchillerNaumann.H"
#include "addToRunTimeSelectionTable.H"

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

namespace Foam
{
    defineTypeNameAndDebug(SchillerNaumann, 0);

    addToRunTimeSelectionTable
    (
        dragModel,
        SchillerNaumann,
        dictionary
    );
}


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

Foam::SchillerNaumann::SchillerNaumann
(
    const dictionary& interfaceDict,
    const volScalarField& alpha,
    const phaseModel& phasea,
    const phaseModel& phaseb
)
:
    dragModel(interfaceDict, alpha, phasea, phaseb)
{}


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

Foam::SchillerNaumann::~SchillerNaumann()
{}


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

Foam::tmp<Foam::volScalarField> Foam::SchillerNaumann::K
(
    const volScalarField& Ur
) const
{
    volScalarField Re = max(Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3));
    
    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();
}


// ************************************************************************* //
SchillerNaumannPatched.C (2,421 bytes)   

albertop

2011-05-25 10:45

reporter  

WenYuPatched.C (2,449 bytes)    
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "WenYu.H"
#include "addToRunTimeSelectionTable.H"

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

namespace Foam
{
    defineTypeNameAndDebug(WenYu, 0);

    addToRunTimeSelectionTable
    (
        dragModel,
        WenYu,
        dictionary
    );
}


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

Foam::WenYu::WenYu
(
    const dictionary& interfaceDict,
    const volScalarField& alpha,
    const phaseModel& phasea,
    const phaseModel& phaseb
)
:
    dragModel(interfaceDict, alpha, phasea, phaseb)
{}


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

Foam::WenYu::~WenYu()
{}


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

Foam::tmp<Foam::volScalarField> Foam::WenYu::K
(
    const volScalarField& Ur
) const
{
    volScalarField beta = max(scalar(1) - alpha_, scalar(1.0e-6));
    volScalarField bp = pow(beta, -2.65);

    volScalarField Re = max(Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3));
    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();
}


// ************************************************************************* //
WenYuPatched.C (2,449 bytes)   

albertop

2011-05-25 10:45

reporter  

SyamlalOBrienPatched.C (2,694 bytes)    
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "SyamlalOBrien.H"
#include "addToRunTimeSelectionTable.H"

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

namespace Foam
{
    defineTypeNameAndDebug(SyamlalOBrien, 0);

    addToRunTimeSelectionTable
    (
        dragModel,
        SyamlalOBrien,
        dictionary
    );
}


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

Foam::SyamlalOBrien::SyamlalOBrien
(
    const dictionary& interfaceDict,
    const volScalarField& alpha,
    const phaseModel& phasea,
    const phaseModel& phaseb
)
:
    dragModel(interfaceDict, alpha, phasea, phaseb)
{}


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

Foam::SyamlalOBrien::~SyamlalOBrien()
{}


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

Foam::tmp<Foam::volScalarField> Foam::SyamlalOBrien::K
(
    const volScalarField& Ur
) const
{
    volScalarField beta = max(scalar(1) - alpha_, scalar(1.0e-6));
    volScalarField A = pow(beta, 4.14);
    volScalarField B = 
	neg(beta - 0.85)*(0.8*pow(beta, 1.28))
      + pos(beta - 0.85)*(pow(beta[celli], 2.65));

    volScalarField Re = max(Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3));

    volScalarField Vr = 0.5*
    (
        A - 0.06*Re + sqrt(sqr(0.06*Re) + 0.12*Re*(2.0*B - A) + sqr(A))
    );

    volScalarField Cds = sqr(0.63 + 4.8*sqrt(Vr/Re));

    return 0.75*Cds*phaseb_.rho()*Ur/(phasea_.d()*sqr(Vr));
}


// ************************************************************************* //
SyamlalOBrienPatched.C (2,694 bytes)   

henry

2011-05-25 11:57

manager   ~0000395

Thanks for the bug-report and the patches, I will merge them in today and push the changes into OpenFOAM-1.7.x

henry

2011-05-25 14:55

manager   ~0000396

Resolved by commit 1e7d927e96835885f4824e750bd182f76424f458

albertop

2011-05-25 21:08

reporter   ~0000404

Thanks for committing the patches. Unfortunately, I introduced a small bug with one of the patches. The GidaspowErgunWenYu.C return statement should read

    return
        pos(beta - 0.8)
       *(0.75*Cds*phaseb_.rho()*Ur*bp/phasea_.d())
      + neg(beta - 0.8)
       *(
           150.0*alpha_*phaseb_.nu()*phaseb_.rho()/(sqr(beta*phasea_.d()))
         + 1.75*phaseb_.rho()*Ur/(beta*phasea_.d())
        );

since K is divided by alpha*beta due to the phase-intensive form of the momentum equation. All other patches are fine.

Sorry for the inconvenience.

henry

2011-05-26 08:27

manager   ~0000409

Resolved by commit 40e7c06e5c65306e78e6a6c2c4fc6867ce250161

Issue History

Date Modified Username Field Change
2011-05-25 00:35 albertop New Issue
2011-05-25 00:35 albertop File Added: bed2.tar.gz
2011-05-25 09:57 albertop File Added: GidaspowErgunWenYu.patch
2011-05-25 09:58 albertop Note Added: 0000392
2011-05-25 10:43 albertop Note Added: 0000393
2011-05-25 10:44 albertop File Added: GidaspowSchillerNaumannPatched.C
2011-05-25 10:44 albertop File Added: SchillerNaumannPatched.C
2011-05-25 10:45 albertop File Added: WenYuPatched.C
2011-05-25 10:45 albertop File Added: SyamlalOBrienPatched.C
2011-05-25 11:57 henry Note Added: 0000395
2011-05-25 14:55 henry Note Added: 0000396
2011-05-25 14:55 henry Status new => resolved
2011-05-25 14:55 henry Resolution open => fixed
2011-05-25 14:55 henry Assigned To => henry
2011-05-25 21:08 albertop Note Added: 0000404
2011-05-25 21:08 albertop Status resolved => feedback
2011-05-25 21:08 albertop Resolution fixed => reopened
2011-05-26 08:27 henry Note Added: 0000409
2011-05-26 08:27 henry Status feedback => resolved
2011-05-26 08:27 henry Resolution reopened => fixed