View Issue Details

IDProjectCategoryView StatusLast Update
0002822OpenFOAMBugpublic2018-01-31 16:45
Reporterrolf.sieber@sms-vt.com Assigned Tohenry  
PriorityhighSeveritymajorReproducibilityalways
Status resolvedResolutionfixed 
PlatformGNU/LinuxOSOtherOS Version(please specify)
Summary0002822: Missinig time dimension in Euler Scheme for moving mesh
DescriptionPossible Error in the file EulerDdtSchme.C in function

tmp<GeometricField<Type, fvPatchField, volMesh>>
EulerDdtScheme<Type>::fvcDdt
(
    const volScalarField& alpha,
    const volScalarField& rho,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)

In case of a moving mesh rDeltaT.value() is used to calculate the time derivative and not rDeltaT like in non moving mesh. This leads to the fact the the dimesion of the returned class is not correct and the program stops with the message

    [h.water[1 -1 -3 0 0 0 0] ] + [ddt(alpha.water,thermo:rho.water,K.water)[1 -1 -2 0 0 0 0] ]

during the calculation of the time derivative of the kinetic energy.
TagsNo tags attached.

Activities

henry

2018-01-31 15:20

manager   ~0009235

Steps to reproduce?

rolf.sieber@sms-vt.com

2018-01-31 15:50

reporter   ~0009236

There is no way in standard branch to reproduce the effect, since no twoPhaseEulerFoam with dynamic mesh exists.

I already implemented dynamic Mesh in twoPhaseEulerFoam in OpenFoam 4.x and got the same error. However I couldn't find the reason.

Now I try to implement dynamic Mesh in reactingTwoPhaseEulerFoam of OpenFoam 5.x and again. Since I'm more familiar with the code I guess I found the reason.

I attached the file AnisothermalPhaseModel.C, where you can find my workaround. (Marked with "RS"). Hope this helps.
AnisothermalPhaseModel.C (5,172 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2015-2017 OpenFOAM Foundation
     \\/     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 "AnisothermalPhaseModel.H"
#include "phaseSystem.H"

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

template<class BasePhaseModel>
Foam::AnisothermalPhaseModel<BasePhaseModel>::AnisothermalPhaseModel
(
    const phaseSystem& fluid,
    const word& phaseName,
    const label index
)
:
    BasePhaseModel(fluid, phaseName, index),
    K_
    (
        IOobject
        (
            IOobject::groupName("K", this->name()),
            fluid.mesh().time().timeName(),
            fluid.mesh()
        ),
        fluid.mesh(),
        dimensionedScalar("K", sqr(dimVelocity), scalar(0))
    )
{}


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

template<class BasePhaseModel>
Foam::AnisothermalPhaseModel<BasePhaseModel>::~AnisothermalPhaseModel()
{}


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

template<class BasePhaseModel>
bool Foam::AnisothermalPhaseModel<BasePhaseModel>::compressible() const
{
    return !this->thermo().incompressible();
}


template<class BasePhaseModel>
void Foam::AnisothermalPhaseModel<BasePhaseModel>::correctKinematics()
{
    BasePhaseModel::correctKinematics();
    K_ = 0.5*magSqr(this->U());
}


template<class BasePhaseModel>
void Foam::AnisothermalPhaseModel<BasePhaseModel>::correctThermo()
{
    BasePhaseModel::correctThermo();

    this->thermo_->correct();
}


template<class BasePhaseModel>
Foam::tmp<Foam::volScalarField>
Foam::AnisothermalPhaseModel<BasePhaseModel>::filterPressureWork
(
    const tmp<volScalarField>& pressureWork
) const
{
    const volScalarField& alpha = *this;

    scalar pressureWorkAlphaLimit =
        this->thermo_->lookupOrDefault("pressureWorkAlphaLimit", 0.0);

    if (pressureWorkAlphaLimit > 0)
    {
        return
        (
            max(alpha - pressureWorkAlphaLimit, scalar(0))
           /max(alpha - pressureWorkAlphaLimit, pressureWorkAlphaLimit)
        )*pressureWork;
    }
    else
    {
        return pressureWork;
    }
}


template<class BasePhaseModel>
Foam::tmp<Foam::fvScalarMatrix>
Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
{
    const volScalarField& alpha = *this;
    const volVectorField& U = this->U();
    const surfaceScalarField& alphaPhi = this->alphaPhi();
    const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi();

    const volScalarField& contErr(this->continuityError());

    const volScalarField alphaEff(this->turbulence().alphaEff());

    volScalarField& he = this->thermo_->he();

    //RS, 31.01.18: Workaround for possibleBug in Eulerscheme
    //              Damit wird die buggy function für die Zeitableitung umgangen
    volScalarField KinEn("KinEn",alpha*this->rho()*K_);
    //RS, 31.01.18: Ende


    tmp<fvScalarMatrix> tEEqn
    (
        fvm::ddt(alpha, this->rho(), he)
      + fvm::div(alphaRhoPhi, he)
      - fvm::Sp(contErr, he)

    //RS, 31.01.18: Workaround for possibleBug in Eulerscheme => Now we are using a different function
    //      + fvc::ddt(alpha, this->rho(), K_) + fvc::div(alphaRhoPhi, K_)
      + fvc::ddt(KinEn) + fvc::div(alphaRhoPhi, K_)
    //RS, 31.01.18: Ende
      - contErr*K_

      - fvm::laplacian
        (
            fvc::interpolate(alpha)
           *fvc::interpolate(alphaEff),
            he
        )
     ==
        alpha*this->Qdot()
    );

    // Add the appropriate pressure-work term
    if (he.name() == this->thermo_->phasePropertyName("e"))
    {
        tEEqn.ref() += filterPressureWork
        (
            fvc::div(fvc::absolute(alphaPhi, alpha, U), this->thermo().p())
          + this->thermo().p()*fvc::ddt(alpha)
        );
    }
    else if (this->thermo_->dpdt())
    {
        tEEqn.ref() -= filterPressureWork(alpha*this->fluid().dpdt());
    }

    return tEEqn;
}


template<class BasePhaseModel>
const Foam::volScalarField&
Foam::AnisothermalPhaseModel<BasePhaseModel>::K() const
{
    return K_;
}


// ************************************************************************* //
AnisothermalPhaseModel.C (5,172 bytes)   

henry

2018-01-31 16:45

manager   ~0009237

Resolved in OpenFOAM-5.x by commit 7bf7ed54c6ed332752b1a8647adc25c6a108f43a

Resolved in OpenFOAM-dev by commit 1a0d663977bd05fb1af23ddbd6813435bc5cb12f

Issue History

Date Modified Username Field Change
2018-01-31 15:05 rolf.sieber@sms-vt.com New Issue
2018-01-31 15:20 henry Note Added: 0009235
2018-01-31 15:50 rolf.sieber@sms-vt.com File Added: AnisothermalPhaseModel.C
2018-01-31 15:50 rolf.sieber@sms-vt.com Note Added: 0009236
2018-01-31 16:45 henry Assigned To => henry
2018-01-31 16:45 henry Status new => resolved
2018-01-31 16:45 henry Resolution open => fixed
2018-01-31 16:45 henry Fixed in Version => 5.x
2018-01-31 16:45 henry Note Added: 0009237