View Issue Details

IDProjectCategoryView StatusLast Update
0002456OpenFOAMBugpublic2020-11-23 12:08
Reporterfsch1 Assigned Tohenry  
PrioritynormalSeveritymajorReproducibilityalways
Status resolvedResolutionfixed 
Fixed in Versiondev 
Summary0002456: curvatureSeparation in 'Hot Boxes'-tutorial doesn't work
DescriptionHello,
in the hotBoxes-tutorial of the reactingParcelFilmFoam the curvatureSeparation doesn't work.

Run the case, e.g 5s. In result-file ./5/uniform/regionModels/wallFilmRegion/wallFilmRegionOutputProperties all the injectedMass is due to drippingInjection and 0 for curvatureSeparation.

In addition, the wallFilm-height deltaf gets to high for my understanding. And then the minimum temperature in the thermoSingleLayer/wallFilm goes down to 270K instead staying at more less 300K.
I think it's a problem for the model if the wallFilm gets to high. But I think it's caused by the water not separating at the bottom of the boxes...
Steps To ReproduceRun hotBoxes-tutorial
TagscurvatureSeparation

Activities

fsch1

2017-02-10 13:08

reporter   ~0007734

I've added a small testcase.

the curvatureSeperation is the only injectionModel activated. And there are some parcels generated. So curvatureSeperation seems to work.

But in 0.05 uniform/regionModels/wallFilmRegion/wallFilmRegionOutputProperties there is no count for the injected parcels by curvatureSeperation...


Run:
Allpre
reactingParcelFilmFoam

Note: the numeric instability is due to the high wallFilm height deltaf in beginning for demonstrating the effect in the first timesteps
bugCurvatureSeperation.zip (3,499,933 bytes)

shildenbrand

2017-02-10 14:58

reporter   ~0007737

I can confirm this problem with curvatureSeparation. I also attached a different testcase which shows this behaviour: When using curvatureSeparation only, the case gets instable and temperatures behave stange.
When using another model (BrunDripping or Dripping) the case stays stable.
I don't think it has something to do with the film height but rather with the model itself.

fsch1

2019-06-21 13:40

reporter   ~0010524

Hello again,
after a long time I took a look into this case again and updated it to work with the newer OpenFoam-dev versions.
In short it's a box with a water film on it, which starts dry by the water running of due to gravity.

To run more stable all the temperatures are the same (293 K) and constrained via fvOptions and the phaseChangeModel is deactivated.

1. Bug
Shildenbrands point is correct regarding the instability of the curvatureSeparation-model. I think it is due to this two lines in the curvatureSeparation.C:
...
 const volScalarField& delta = film.delta();
...
diameterToInject = separated*delta;
...

As far as I understand the code, if droplets are separating, the diameter of these droplets are the same as the film height.
So this does not make sense.
Additionally if the film height is really low, the 1e8 to 1e11 particle are added into one single parcel. This is also the point when the simulation becomes instable (in the attached case shortly after 0.07s).

I think it's necessary to implement some checks (like in drippingInjection.C, starting line 129), that droplets only can separate, if a minMass is overcome. For the diameter I don't know what's the best solution.


2.Bug
The mass injected by the curvatureSeparation-Model continues to be neglected in the outputProperties.
As you can see in the wallFilmRegionOutputProperties below injectedMass is zero, even more ~1500 parcels detached the film.

...
    location "0.07/uniform/regionModels/wallFilmRegion";
    object wallFilmRegionOutputProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

addedMassTotal 0.006046627137;

injectionModel
{
    curvatureSeparation
    {
        injectedMass 0;
    }
}



Hope this helps in locating these bugs.

fsch1

2019-07-01 13:13

reporter   ~0010537

Regarding 1.)
I have attached a proposed patch of the curvatureSeparation.
The patch introduces a stable film thickness (similar as in the two other existing injectionModels brunDrippingInjection/ drippingInjection) which have to be exceeded for the model to activate the separation.

With a film tickness deltaStable = 0.0001 the attached test case as well as the Hot Boxes - Tutorial stay numerically stable.
curvatureSeparation.C (10,469 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Copyright (C) 2011-2018 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 "curvatureSeparation.H"
#include "addToRunTimeSelectionTable.H"
#include "fvMesh.H"
#include "Time.H"
#include "volFields.H"
#include "kinematicSingleLayer.H"
#include "surfaceInterpolate.H"
#include "fvcDiv.H"
#include "fvcGrad.H"
#include "stringListOps.H"
#include "cyclicPolyPatch.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{
namespace regionModels
{
namespace surfaceFilmModels
{

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

defineTypeNameAndDebug(curvatureSeparation, 0);
addToRunTimeSelectionTable
(
    injectionModel,
    curvatureSeparation,
    dictionary
);

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

tmp<volScalarField> curvatureSeparation::calcInvR1
(
    const volVectorField& U
) const
{
    // method 1
/*
    tmp<volScalarField> tinvR1
    (
        volScalarField::New("invR1", fvc::div(film().nHat()))
    );
*/

    // method 2
    dimensionedScalar smallU("smallU", dimVelocity, rootVSmall);
    volVectorField UHat(U/(mag(U) + smallU));
    tmp<volScalarField> tinvR1
    (
        volScalarField::New("invR1", UHat & (UHat & gradNHat_))
    );


    scalarField& invR1 = tinvR1.ref().primitiveFieldRef();

    // apply defined patch radii
    const scalar rMin = 1e-6;
    const fvMesh& mesh = film().regionMesh();
    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
    forAll(definedPatchRadii_, i)
    {
        label patchi = definedPatchRadii_[i].first();
        scalar definedInvR1 = 1.0/max(rMin, definedPatchRadii_[i].second());
        UIndirectList<scalar>(invR1, pbm[patchi].faceCells()) = definedInvR1;
    }

    // filter out large radii
    const scalar rMax = 1e6;
    forAll(invR1, i)
    {
        if (mag(invR1[i]) < 1/rMax)
        {
            invR1[i] = -1.0;
        }
    }

    if (debug && mesh.time().writeTime())
    {
        tinvR1().write();
    }

    return tinvR1;
}


tmp<scalarField> curvatureSeparation::calcCosAngle
(
    const surfaceScalarField& phi
) const
{
    const fvMesh& mesh = film().regionMesh();
    const vectorField nf(mesh.Sf()/mesh.magSf());
    const unallocLabelList& own = mesh.owner();
    const unallocLabelList& nbr = mesh.neighbour();

    scalarField phiMax(mesh.nCells(), -great);
    scalarField cosAngle(mesh.nCells(), 0.0);
    forAll(nbr, facei)
    {
        label cellO = own[facei];
        label cellN = nbr[facei];

        if (phi[facei] > phiMax[cellO])
        {
            phiMax[cellO] = phi[facei];
            cosAngle[cellO] = -gHat_ & nf[facei];
        }
        if (-phi[facei] > phiMax[cellN])
        {
            phiMax[cellN] = -phi[facei];
            cosAngle[cellN] = -gHat_ & -nf[facei];
        }
    }

    forAll(phi.boundaryField(), patchi)
    {
        const fvsPatchScalarField& phip = phi.boundaryField()[patchi];
        const fvPatch& pp = phip.patch();
        const labelList& faceCells = pp.faceCells();
        const vectorField nf(pp.nf());
        forAll(phip, i)
        {
            label celli = faceCells[i];
            if (phip[i] > phiMax[celli])
            {
                phiMax[celli] = phip[i];
                cosAngle[celli] = -gHat_ & nf[i];
            }
        }
    }
/*
    // correction for cyclics - use cyclic pairs' face normal instead of
    // local face normal
    const fvBoundaryMesh& pbm = mesh.boundary();
    forAll(phi.boundaryField(), patchi)
    {
        if (isA<cyclicPolyPatch>(pbm[patchi]))
        {
            const scalarField& phip = phi.boundaryField()[patchi];
            const vectorField nf(pbm[patchi].nf());
            const labelList& faceCells = pbm[patchi].faceCells();
            const label sizeBy2 = pbm[patchi].size()/2;

            for (label face0=0; face0<sizeBy2; face0++)
            {
                label face1 = face0 + sizeBy2;
                label cell0 = faceCells[face0];
                label cell1 = faceCells[face1];

                // flux leaving half 0, entering half 1
                if (phip[face0] > phiMax[cell0])
                {
                    phiMax[cell0] = phip[face0];
                    cosAngle[cell0] = -gHat_ & -nf[face1];
                }

                // flux leaving half 1, entering half 0
                if (-phip[face1] > phiMax[cell1])
                {
                    phiMax[cell1] = -phip[face1];
                    cosAngle[cell1] = -gHat_ & nf[face0];
                }
            }
        }
    }
*/
    // checks
    if (debug && mesh.time().writeTime())
    {
        volScalarField volCosAngle
        (
            IOobject
            (
                "cosAngle",
                mesh.time().timeName(),
                mesh,
                IOobject::NO_READ
            ),
            mesh,
            dimensionedScalar(dimless, 0),
            zeroGradientFvPatchScalarField::typeName
        );
        volCosAngle.primitiveFieldRef() = cosAngle;
        volCosAngle.correctBoundaryConditions();
        volCosAngle.write();
    }

    return max(min(cosAngle, scalar(1)), scalar(-1));
}


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

curvatureSeparation::curvatureSeparation
(
    surfaceFilmRegionModel& film,
    const dictionary& dict
)
:
    injectionModel(type(), film, dict),
    gradNHat_(fvc::grad(film.nHat())),
    deltaByR1Min_(coeffDict_.lookupOrDefault<scalar>("deltaByR1Min", 0.0)),
    definedPatchRadii_(),
    magG_(mag(film.g().value())),
    gHat_(Zero),
    deltaStable_(coeffDict_.lookupOrDefault("deltaStable", scalar(0)))
{
    if (magG_ < rootVSmall)
    {
        FatalErrorInFunction
            << "Acceleration due to gravity must be non-zero"
            << exit(FatalError);
    }

    gHat_ = film.g().value()/magG_;

    List<Tuple2<word, scalar>> prIn(coeffDict_.lookup("definedPatchRadii"));
    const wordList& allPatchNames = film.regionMesh().boundaryMesh().names();

    DynamicList<Tuple2<label, scalar>> prData(allPatchNames.size());

    labelHashSet uniquePatchIDs;

    forAllReverse(prIn, i)
    {
        labelList patchIDs = findStrings(prIn[i].first(), allPatchNames);
        forAll(patchIDs, j)
        {
            const label patchi = patchIDs[j];

            if (!uniquePatchIDs.found(patchi))
            {
                const scalar radius = prIn[i].second();
                prData.append(Tuple2<label, scalar>(patchi, radius));

                uniquePatchIDs.insert(patchi);
            }
        }
    }

    definedPatchRadii_.transfer(prData);
}


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

curvatureSeparation::~curvatureSeparation()
{}


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

void curvatureSeparation::correct
(
    scalarField& availableMass,
    scalarField& massToInject,
    scalarField& diameterToInject
)
{
    const kinematicSingleLayer& film =
        refCast<const kinematicSingleLayer>(this->film());
    const fvMesh& mesh = film.regionMesh();

    const volScalarField& delta = film.delta();
    const volVectorField& U = film.U();
    const surfaceScalarField& phi = film.phi();
    const volScalarField& rho = film.rho();
    const scalarField magSqrU(magSqr(film.U()));
    const volScalarField& sigma = film.sigma();

    const scalarField invR1(calcInvR1(U));
    const scalarField cosAngle(calcCosAngle(phi));

    // calculate force balance
    const scalar Fthreshold = 1e-10;
    scalarField Fnet(mesh.nCells(), 0.0);
    scalarField separated(mesh.nCells(), 0.0);
    forAll(invR1, i)
    {
        if ((invR1[i] > 0) && (delta[i]*invR1[i] > deltaByR1Min_))
        {
            scalar R1 = 1.0/(invR1[i] + rootVSmall);
            scalar R2 = R1 + delta[i];

            // inertial force
            scalar Fi = -delta[i]*rho[i]*magSqrU[i]*72.0/60.0*invR1[i];

            // body force
            scalar Fb =
              - 0.5*rho[i]*magG_*invR1[i]*(sqr(R1) - sqr(R2))*cosAngle[i];

            // surface force
            scalar Fs = sigma[i]/R2;

            Fnet[i] = Fi + Fb + Fs;
 
            if ((Fnet[i] + Fthreshold < 0) && (delta[i] > deltaStable_))
            {
		separated[i] = 1.0;		
            }
        }
    }

    // inject all available mass
    massToInject = separated*availableMass;
    diameterToInject = separated*delta;
    availableMass -= separated*availableMass;

    addToInjectedMass(sum(separated*availableMass));

    if (debug && mesh.time().writeTime())
    {
        volScalarField volFnet
        (
            IOobject
            (
                "Fnet",
                mesh.time().timeName(),
                mesh,
                IOobject::NO_READ
            ),
            mesh,
            dimensionedScalar(dimForce, 0),
            zeroGradientFvPatchScalarField::typeName
        );
        volFnet.primitiveFieldRef() = Fnet;
        volFnet.correctBoundaryConditions();
        volFnet.write();
    }

    injectionModel::correct();
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace surfaceFilmModels
} // End namespace regionModels
} // End namespace Foam

// ************************************************************************* //
curvatureSeparation.C (10,469 bytes)   
curvatureSeparation.H (4,520 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Copyright (C) 2011-2019 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/>.

Class
    Foam::regionModels::surfaceFilmModels::curvatureSeparation

Description
    Curvature film separation model

    Assesses film curvature via the mesh geometry and calculates a force
    balance of the form:

        F_sum = F_inertial + F_body + F_surface

    If F_sum < 0, the film separates. Similarly, if F_sum > 0 the film will
    remain attached.

    Based on description given by
        Owen and D. J. Ryley. The flow of thin liquid films around corners.
        International Journal of Multiphase Flow, 11(1):51-62, 1985.


SourceFiles
    curvatureSeparation.C

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

#ifndef curvatureSeparation_H
#define curvatureSeparation_H

#include "injectionModel.H"
#include "surfaceFields.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{
namespace regionModels
{
namespace surfaceFilmModels
{

/*---------------------------------------------------------------------------*\
                   Class curvatureSeparation Declaration
\*---------------------------------------------------------------------------*/

class curvatureSeparation
:
    public injectionModel
{
protected:

    // Protected data

        //- Gradient of surface normals
        volTensorField gradNHat_;

        //- Minimum gravity driven film thickness (non-dimensionalised delta/R1)
        scalar deltaByR1Min_;

        //- List of radii for patches - if patch not defined, radius
        // calculated based on mesh geometry
        List<Tuple2<label, scalar>> definedPatchRadii_;

        //- Magnitude of gravity vector
        scalar magG_;

        //- Direction of gravity vector
        vector gHat_;

	//- Stable film thickness - separation only possible if thickness
	//  exceeds this threshold value
	scalar deltaStable_;


    // Protected Member Functions

        //- Calculate local (inverse) radius of curvature
        tmp<volScalarField> calcInvR1(const volVectorField& U) const;

        //- Calculate the cosine of the angle between gravity vector and
        //  cell out flow direction
        tmp<scalarField> calcCosAngle(const surfaceScalarField& phi) const;


public:

    //- Runtime type information
    TypeName("curvatureSeparation");


    // Constructors

        //- Construct from surface film model
        curvatureSeparation
        (
            surfaceFilmRegionModel& film,
            const dictionary& dict
        );

        //- Disallow default bitwise copy construction
        curvatureSeparation(const curvatureSeparation&) = delete;


    //- Destructor
    virtual ~curvatureSeparation();


    // Member Functions

        // Evolution

            //- Correct
            virtual void correct
            (
                scalarField& availableMass,
                scalarField& massToInject,
                scalarField& diameterToInject
            );


    // Member Operators

        //- Disallow default bitwise assignment
        void operator=(const curvatureSeparation&) = delete;
};


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace surfaceFilmModels
} // End namespace regionModels
} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //
curvatureSeparation.H (4,520 bytes)   

henry

2020-11-23 12:08

manager   ~0011738

Resolved by commit d7d1221cd47833c0e12fbbe5e38bc5f3db95a467

Issue History

Date Modified Username Field Change
2017-02-09 13:43 fsch1 New Issue
2017-02-09 13:43 fsch1 Tag Attached: curvatureSeparation
2017-02-10 13:08 fsch1 File Added: bugCurvatureSeperation.zip
2017-02-10 13:08 fsch1 Note Added: 0007734
2017-02-10 14:58 shildenbrand File Added: testcase_curvatureSeparation.tgz
2017-02-10 14:58 shildenbrand Note Added: 0007737
2019-06-21 13:40 fsch1 File Added: V03_TestCase_CurvatureSeparation.zip
2019-06-21 13:40 fsch1 Note Added: 0010524
2019-07-01 13:13 fsch1 File Added: curvatureSeparation.C
2019-07-01 13:13 fsch1 File Added: curvatureSeparation.H
2019-07-01 13:13 fsch1 Note Added: 0010537
2020-11-23 12:08 henry Assigned To => henry
2020-11-23 12:08 henry Status new => resolved
2020-11-23 12:08 henry Resolution open => fixed
2020-11-23 12:08 henry Fixed in Version => dev
2020-11-23 12:08 henry Note Added: 0011738