View Issue Details

IDProjectCategoryView StatusLast Update
0002713OpenFOAMPatchpublic2017-11-07 13:56
Reporterbjoern.pfeiffelmann Assigned Tohenry  
PrioritynormalSeveritymajorReproducibilityalways
Status resolvedResolutionfixed 
PlatformGNU/LinuxOSOpenSuSEOS Version42.3
Fixed in Versiondev 
Summary0002713: Harmonic surface interpolation scheme predict the boundary face values wrong
DescriptionAs I have posted already here: https://www.cfd-online.com/Forums/openfoam-solving/192948-chtmultiregionsimplefoam-nonlinear-heat-conduction-wrong-temperature-prediction.html
is the number of cells required to predict a nonlinear diffusion problem with the chtMultiRegionSimpleFoam solver quite high.

To improve this behavior the harmonic surface interpolation scheme for interpolating the diffusion coefficient seems to be the right choice as described by e.g. “Numerical heat transfer and fluid flow. S. V. Patankar (1980)”
But there is a bug in the interpolation of the boundary patch values in the harmonic scheme, more precisely in the reverse Linear scheme. As can been seen in the figure attached, I would expect that the interpolated boundary field value is equal to the internal field value and not equal to the boundary field if the harmonic scheme is applied.


We have created a patch for the problem which is attached including the validation case. We modified the reverse Linear scheme, so that it corrects the interpolated boundary field afterwards. This helped to be able to predict the solution only with 3 cells but is rather less elegant.

What is the reason to ignore the weighting factor in surfaceInterpolationScheme.C:308 for the boundary patch? Shouldn’t be the weighting factor dependent on the interpolation scheme instead of setting it equal to 1 like in fvPatch.C:152

best regards
Björn
Steps To Reproducedownload harmonicPatch.zip
unzip
run Allrun
Tagsharmonic, Heat transfer, Surface Interpolation

Activities

bjoern.pfeiffelmann

2017-10-06 13:35

reporter  

harmonicPatch.zip (315,144 bytes)

bjoern.pfeiffelmann

2017-10-06 13:36

reporter  

harmonicPatch.jpg (77,410 bytes)   
harmonicPatch.jpg (77,410 bytes)   

bjoern.pfeiffelmann

2017-10-06 13:36

reporter  

OF_vs_ANALYTIC_T.png (6,929 bytes)   
OF_vs_ANALYTIC_T.png (6,929 bytes)   

bjoern.pfeiffelmann

2017-10-06 13:37

reporter  

henry

2017-10-06 15:04

manager   ~0008825

Adapting the harmonic interpolation to override the boundary value is not appropriate in general; for example if it were used for a transported field the value on the boundary is the correct value for convective transport into the domain.

What you really need is a wall-function BC for the diffusion coefficient which sets the value on the boundary appropriate for the transport between the boundary and cell-centre.

StephanG

2017-10-06 18:21

reporter   ~0008829

I have cleaned up your case and added an automatic graph creation to it. Might make resolving this issue faster. @henry wouldn't it make sense to add such an alpha correction to every (non coupled) boundary by default?
wall1D.tar.gz (3,950 bytes)

henry

2017-10-06 18:31

manager   ~0008830

What you really need is a wall-function BC for the diffusion coefficient which sets the value on the boundary appropriate for the transport between the boundary and cell-centre.

StephanG

2017-10-06 19:03

reporter   ~0008831

Yes, but what I meant is: Shouldn't this be part of thermo.correct()? Or am i missing something?

henry

2017-10-06 20:12

manager   ~0008832

It would be a boundary condition and updated as a boundary condition.

bjoern.pfeiffelmann

2017-10-09 07:53

reporter   ~0008833

I think that it is not necessary to develop an additional boundary condition. The problem is that the harmonic scheme is not implemented consistently (in the internal field yes, but on the boundary not). The present scheme aim to correct this. I think that this does not have any implications for the convective transport, since the harmonic scheme is not used for the convective therm anyway.
For my understanding linear and reverseLinear are opposite in the sense that the linear scheme weights the cell center higher which is closer to the face center for which the field is interpolated. reverseLinear weights the cell center higher which is further away from the face center. But why should this be only valid for internal fields or coupled boundaries? Since the boundary patch field (not interpolated) is positioned exactly on the boundary face the boundary patch field shouldn’t influence the interpolated boundary patch field in the case of using the reverseLinear scheme. So the result is that the interpolated boundary field is equal to the internal field because in comparison with the boundary patch field it has a distance to the boundary face.

Here is my proposal to implement this behavior
https://github.com/OpenFOAM/OpenFOAM-dev/pull/13

Implemented in that way it should not influence any other scheme, because the weighting for the boundary for all other schemes is 1. It also includes additional flexibility in the implementation of other schemes.

henry

2017-10-09 08:20

manager   ~0008834

Last edited: 2017-10-09 08:51

The harmonic interpolation is run-time selectable and could be used for convective transport or anything else and should not override BCs which are fixed or specifically evaluated in some way. Boundary values are boundary values and not interpolated unless they are coupled with another domain, e.g. processor boundaries which which case they are not real boundaries anyway.

If you want to extrapolate the internal field values to the boundary for whatever reason you can simply apply the 'zeroGradient' BC, you do not need to change the interpolation scheme into an extrapolation scheme to do this.

bjoern.pfeiffelmann

2017-10-09 09:48

reporter   ~0008835

>The harmonic interpolation is run-time selectable and could be used for convective transport or anything else and should not override BCs which are fixed or specifically evaluated in some way. Boundary values are boundary values and not interpolated unless they are coupled with another domain, e.g. processor boundaries which which case they are not real boundaries anyway.

I understand the doubts, but there is also a disadvantage in creating separate boundary conditions for the harmonic interpolation scheme: When the user on run time level specifies the harmonic interpolation scheme, without the (new) harmonic boundary condition. In that case the intended improvement would be ineffective.

>If you want to extrapolate the internal field values to the boundary for whatever reason you can simply apply the 'zeroGradient' BC, you do not need to change the interpolation scheme into an extrapolation scheme to do this.

My aim, and I think that it is usually the aim when using the harmonic interpolation, is the surface interpolation of the diffusion coefficient but not the interpolation of the field variable. Is it possible, at all, to apply 'zeroGradient' BC for a diffusion coefficient?

henry

2017-10-09 10:04

manager   ~0008836

> Is it possible, at all, to apply 'zeroGradient' BC for a diffusion coefficient?

I have not checked this but if it is not currently possible it should be made possible rather then compromising the design of the interpolation schemes.

StephanG

2017-10-09 14:23

reporter   ~0008843

Currently the only option for solids is heSolidThermo. Which uses kappa (NO_READ NO_WRITE) to calculate alpha which the solver uses (from basic thermo).thermo:alpha is again set to NO_READ and NO_WRITE. The thermo.correct() or calculate() function updates the values for the temperature dependence. Hence my statement above. From my perspective the change should still be made there.

henry

2017-10-09 14:28

manager   ~0008844

Agreed the interpolation scheme should not manipulate boundary values, this is what boundary conditions are for, the thermodynamics should be generalized to allow specialized BCs to be applied to kappa.

bjoern.pfeiffelmann

2017-10-10 07:33

reporter   ~0008846

When implemented as a specialized BC, I would prefer a solution which is not limited to heat transfer problems or a wall boundary condition. This issue should occur whenever a diffusion flux is approximated at the boundary with varying diffusion coefficients e.g. species concentration at a velocity inlet.

henry

2017-10-10 08:06

manager   ~0008847

Can you provide a patch which adds run-time selected BCs to kappa and diffusivity?

bjoern.pfeiffelmann

2017-10-10 13:46

reporter   ~0008851

I will try it, but it will take some time

henry

2017-10-13 09:31

manager   ~0008865

Please reopen when you can have completed the patch which adds run-time selected BCs to kappa and diffusivity.

Thanks

bjoern.pfeiffelmann

2017-10-20 12:27

reporter   ~0008907

Here my patch for solid thermo:alpha:

https://github.com/OpenFOAM/OpenFOAM-dev/pull/15

and the modified test case attached.

Unfortunately, same problem will occur when the harmonic interpolation scheme is applied to any other diffusion coefficient.

bjoern.pfeiffelmann

2017-10-20 12:28

reporter  

wall1Dharmonic.tar.gz (3,829 bytes)

henry

2017-10-20 12:52

manager   ~0008908

> Unfortunately, same problem will occur when the harmonic interpolation scheme is applied to any other diffusion coefficient.

Sure, have you tried applying the same change to any other diffusion coefficient?

bjoern.pfeiffelmann

2017-10-20 15:23

reporter   ~0008913

No, not jet. Maybe a note can be putted in the description of the harmonic scheme

henry

2017-10-20 15:27

manager   ~0008914

> Maybe a note can be putted in the description of the harmonic scheme

This has nothing to do with the harmonic scheme which only handles interpolation; what you want to do is change the boundary conditions.

bjoern.pfeiffelmann

2017-10-23 07:58

reporter   ~0008919

Are the selection of the interpolation scheme for the diffusion coefficient and the specification of the diffusion coefficient at the boundaries really not related to each other? What about the provided derivation in harmonicPatch.jpg? Are there any wrong assumptions inside? The derivation is based on the description of the books written by Patankar S.V.:

Numerical Heat Transfer and Fluid Flow (1980)
harmonic scheme: pp.44-47, treatment of the boundary cells: pp.68-71

Computation of Conduction and Duct Flow Heat Transfer (1991)
harmonic scheme: pp.36-38, treatment of the boundary cells: pp.43-45, Usage of internal diffusion coefficient as boundary diffusion coefficient: Eq. 5.24 p.77

I thought that you prefer the consistent separation of interpolation and bc (There are good reasons for that), although they are connected to each other. Therefore I provided the last patch with this separation in mind.

bjoern.pfeiffelmann

2017-10-23 08:04

reporter  

basicThermo.C (12,231 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-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 "basicThermo.H"
#include "zeroGradientFvPatchFields.H"
#include "fixedEnergyFvPatchScalarField.H"
#include "gradientEnergyFvPatchScalarField.H"
#include "mixedEnergyFvPatchScalarField.H"
#include "fixedJumpFvPatchFields.H"
#include "fixedJumpAMIFvPatchFields.H"
#include "energyJumpFvPatchScalarField.H"
#include "energyJumpAMIFvPatchScalarField.H"


/* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */

namespace Foam
{
    defineTypeNameAndDebug(basicThermo, 0);
    defineRunTimeSelectionTable(basicThermo, fvMesh);
}

const Foam::word Foam::basicThermo::dictName("thermophysicalProperties");


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

Foam::wordList Foam::basicThermo::heBoundaryBaseTypes()
{
    const volScalarField::Boundary& tbf =
        this->T_.boundaryField();

    wordList hbt(tbf.size(), word::null);

    forAll(tbf, patchi)
    {
        if (isA<fixedJumpFvPatchScalarField>(tbf[patchi]))
        {
            const fixedJumpFvPatchScalarField& pf =
                dynamic_cast<const fixedJumpFvPatchScalarField&>(tbf[patchi]);

            hbt[patchi] = pf.interfaceFieldType();
        }
        else if (isA<fixedJumpAMIFvPatchScalarField>(tbf[patchi]))
        {
            const fixedJumpAMIFvPatchScalarField& pf =
                dynamic_cast<const fixedJumpAMIFvPatchScalarField&>
                (
                    tbf[patchi]
                );

            hbt[patchi] = pf.interfaceFieldType();
        }
    }

    return hbt;
}


Foam::wordList Foam::basicThermo::heBoundaryTypes()
{
    const volScalarField::Boundary& tbf =
        this->T_.boundaryField();

    wordList hbt = tbf.types();

    forAll(tbf, patchi)
    {
        if (isA<fixedValueFvPatchScalarField>(tbf[patchi]))
        {
            hbt[patchi] = fixedEnergyFvPatchScalarField::typeName;
        }
        else if
        (
            isA<zeroGradientFvPatchScalarField>(tbf[patchi])
         || isA<fixedGradientFvPatchScalarField>(tbf[patchi])
        )
        {
            hbt[patchi] = gradientEnergyFvPatchScalarField::typeName;
        }
        else if (isA<mixedFvPatchScalarField>(tbf[patchi]))
        {
            hbt[patchi] = mixedEnergyFvPatchScalarField::typeName;
        }
        else if (isA<fixedJumpFvPatchScalarField>(tbf[patchi]))
        {
            hbt[patchi] = energyJumpFvPatchScalarField::typeName;
        }
        else if (isA<fixedJumpAMIFvPatchScalarField>(tbf[patchi]))
        {
            hbt[patchi] = energyJumpAMIFvPatchScalarField::typeName;
        }
        else if (tbf[patchi].type() == "energyRegionCoupledFvPatchScalarField")
        {
            hbt[patchi] = "energyRegionCoupledFvPatchScalarField";
        }
    }

    return hbt;
}


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

Foam::volScalarField& Foam::basicThermo::lookupOrConstruct
(
    const fvMesh& mesh,
    const char* name
) const
{
    if (!mesh.objectRegistry::foundObject<volScalarField>(name))
    {
        volScalarField* fPtr
        (
            new volScalarField
            (
                IOobject
                (
                    name,
                    mesh.time().timeName(),
                    mesh,
                    IOobject::MUST_READ,
                    IOobject::AUTO_WRITE
                ),
                mesh
            )
        );

        // Transfer ownership of this object to the objectRegistry
        fPtr->store(fPtr);
    }

    return mesh.objectRegistry::lookupObjectRef<volScalarField>(name);
}


Foam::basicThermo::basicThermo
(
    const fvMesh& mesh,
    const word& phaseName
)
:
    IOdictionary
    (
        IOobject
        (
            phasePropertyName(dictName, phaseName),
            mesh.time().constant(),
            mesh,
            IOobject::MUST_READ_IF_MODIFIED,
            IOobject::NO_WRITE
        )
    ),

    phaseName_(phaseName),

    p_(lookupOrConstruct(mesh, "p")),

    T_
    (
        IOobject
        (
            phasePropertyName("T"),
            mesh.time().timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    ),

    alpha_
    (
        IOobject
        (
            phasePropertyName("thermo:alpha"),
            mesh.time().timeName(),
            mesh,
            IOobject::READ_IF_PRESENT,
            IOobject::NO_WRITE
        ),
        mesh,
        dimensionedScalar
        (
            "zero",
            dimensionSet(1, -1, -1, 0, 0),
            Zero
        )
    ),

    dpdt_(lookupOrDefault<Switch>("dpdt", true))
{}


Foam::basicThermo::basicThermo
(
    const fvMesh& mesh,
    const dictionary& dict,
    const word& phaseName
)
:
    IOdictionary
    (
        IOobject
        (
            phasePropertyName(dictName, phaseName),
            mesh.time().constant(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        dict
    ),

    phaseName_(phaseName),

    p_(lookupOrConstruct(mesh, "p")),

    T_
    (
        IOobject
        (
            phasePropertyName("T"),
            mesh.time().timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    ),

    alpha_
    (
        IOobject
        (
            phasePropertyName("thermo:alpha"),
            mesh.time().timeName(),
            mesh,
            IOobject::READ_IF_PRESENT,
            IOobject::NO_WRITE
        ),
        mesh,
        dimensionedScalar
        (
            "zero",
            dimensionSet(1, -1, -1, 0, 0),
            Zero
        )
    )
{}


// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //

Foam::autoPtr<Foam::basicThermo> Foam::basicThermo::New
(
    const fvMesh& mesh,
    const word& phaseName
)
{
    return New<basicThermo>(mesh, phaseName);
}


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

Foam::basicThermo::~basicThermo()
{}


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

const Foam::basicThermo& Foam::basicThermo::lookupThermo
(
    const fvPatchScalarField& pf
)
{
    if (pf.db().foundObject<basicThermo>(dictName))
    {
        return pf.db().lookupObject<basicThermo>(dictName);
    }
    else
    {
        HashTable<const basicThermo*> thermos =
            pf.db().lookupClass<basicThermo>();

        for
        (
            HashTable<const basicThermo*>::iterator iter = thermos.begin();
            iter != thermos.end();
            ++iter
        )
        {
            if
            (
                &(iter()->he().internalField())
              == &(pf.internalField())
            )
            {
                return *iter();
            }
        }
    }

    return pf.db().lookupObject<basicThermo>(dictName);
}


void Foam::basicThermo::validate
(
    const string& app,
    const word& a
) const
{
    if (!(he().name() == phasePropertyName(a)))
    {
        FatalErrorInFunction
            << "Supported energy type is " << phasePropertyName(a)
            << ", thermodynamics package provides " << he().name()
            << exit(FatalError);
    }
}

void Foam::basicThermo::validate
(
    const string& app,
    const word& a,
    const word& b
) const
{
    if
    (
       !(
            he().name() == phasePropertyName(a)
         || he().name() == phasePropertyName(b)
        )
    )
    {
        FatalErrorInFunction
            << "Supported energy types are " << phasePropertyName(a)
            << " and " << phasePropertyName(b)
            << ", thermodynamics package provides " << he().name()
            << exit(FatalError);
    }
}

void Foam::basicThermo::validate
(
    const string& app,
    const word& a,
    const word& b,
    const word& c
) const
{
    if
    (
       !(
            he().name() == phasePropertyName(a)
         || he().name() == phasePropertyName(b)
         || he().name() == phasePropertyName(c)
        )
    )
    {
        FatalErrorInFunction
            << "Supported energy types are " << phasePropertyName(a)
            << ", " << phasePropertyName(b)
            << " and " << phasePropertyName(c)
            << ", thermodynamics package provides " << he().name()
            << exit(FatalError);
    }
}

void Foam::basicThermo::validate
(
    const string& app,
    const word& a,
    const word& b,
    const word& c,
    const word& d
) const
{
    if
    (
       !(
            he().name() == phasePropertyName(a)
         || he().name() == phasePropertyName(b)
         || he().name() == phasePropertyName(c)
         || he().name() == phasePropertyName(d)
        )
    )
    {
        FatalErrorInFunction
            << "Supported energy types are " << phasePropertyName(a)
            << ", " << phasePropertyName(b)
            << ", " << phasePropertyName(c)
            << " and " << phasePropertyName(d)
            << ", thermodynamics package provides " << he().name()
            << exit(FatalError);
    }
}


Foam::wordList Foam::basicThermo::splitThermoName
(
    const word& thermoName,
    const int nCmpt
)
{
    wordList cmpts(nCmpt);

    string::size_type beg=0, end=0, endb=0, endc=0;
    int i = 0;

    while
    (
        (endb = thermoName.find('<', beg)) != string::npos
     || (endc = thermoName.find(',', beg)) != string::npos
    )
    {
        if (endb == string::npos)
        {
            end = endc;
        }
        else if ((endc = thermoName.find(',', beg)) != string::npos)
        {
            end = min(endb, endc);
        }
        else
        {
            end = endb;
        }

        if (beg < end)
        {
            cmpts[i] = thermoName.substr(beg, end-beg);
            cmpts[i++].replaceAll(">","");

            // If the number of number of components in the name
            // is greater than nCmpt return an empty list
            if (i == nCmpt)
            {
                return wordList::null();
            }
        }
        beg = end + 1;
    }

    // If the number of number of components in the name is not equal to nCmpt
    // return an empty list
    if (i + 1 != nCmpt)
    {
        return wordList::null();
    }

    if (beg < thermoName.size())
    {
        cmpts[i] = thermoName.substr(beg, string::npos);
        cmpts[i].replaceAll(">","");
    }

    return cmpts;
}


Foam::volScalarField& Foam::basicThermo::p()
{
    return p_;
}


const Foam::volScalarField& Foam::basicThermo::p() const
{
    return p_;
}


const Foam::volScalarField& Foam::basicThermo::T() const
{
    return T_;
}


Foam::volScalarField& Foam::basicThermo::T()
{
    return T_;
}


const Foam::volScalarField& Foam::basicThermo::alpha() const
{
    return alpha_;
}


const Foam::scalarField& Foam::basicThermo::alpha(const label patchi) const
{
    return alpha_.boundaryField()[patchi];
}


bool Foam::basicThermo::read()
{
    return regIOobject::read();
}


// ************************************************************************* //
basicThermo.C (12,231 bytes)   

bjoern.pfeiffelmann

2017-10-23 08:04

reporter  

heSolidThermo.C (8,236 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2016 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 "heSolidThermo.H"
#include "volFields.H"

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

template<class BasicSolidThermo, class MixtureType>
void Foam::heSolidThermo<BasicSolidThermo, MixtureType>::calculate()
{
    scalarField& TCells = this->T_.primitiveFieldRef();

    const scalarField& hCells = this->he_;
    const scalarField& pCells = this->p_;
    scalarField& rhoCells = this->rho_.primitiveFieldRef();
    scalarField& alphaCells = this->alpha_.primitiveFieldRef();

    forAll(TCells, celli)
    {
        const typename MixtureType::thermoType& mixture_ =
            this->cellMixture(celli);

        const typename MixtureType::thermoType& volMixture_ =
            this->cellVolMixture(pCells[celli], TCells[celli], celli);

        TCells[celli] = mixture_.THE
        (
            hCells[celli],
            pCells[celli],
            TCells[celli]
        );

        rhoCells[celli] = volMixture_.rho(pCells[celli], TCells[celli]);

        alphaCells[celli] =
            volMixture_.kappa(pCells[celli], TCells[celli])
            /
            mixture_.Cpv(pCells[celli], TCells[celli]);
    }

    volScalarField::Boundary& pBf =
        this->p_.boundaryFieldRef();

    volScalarField::Boundary& TBf =
        this->T_.boundaryFieldRef();

    volScalarField::Boundary& rhoBf =
        this->rho_.boundaryFieldRef();

    volScalarField::Boundary& heBf =
        this->he().boundaryFieldRef();

    volScalarField::Boundary& alphaBf =
        this->alpha_.boundaryFieldRef();

    forAll(this->T_.boundaryField(), patchi)
    {
        fvPatchScalarField& pp = pBf[patchi];
        fvPatchScalarField& pT = TBf[patchi];
        fvPatchScalarField& prho = rhoBf[patchi];
        fvPatchScalarField& phe = heBf[patchi];
        fvPatchScalarField& palpha = alphaBf[patchi];

        if (pT.fixesValue())
        {
            forAll(pT, facei)
            {
                const typename MixtureType::thermoType& mixture_ =
                    this->patchFaceMixture(patchi, facei);

                const typename MixtureType::thermoType& volMixture_ =
                    this->patchFaceVolMixture
                    (
                        pp[facei],
                        pT[facei],
                        patchi,
                        facei
                    );


                phe[facei] = mixture_.HE(pp[facei], pT[facei]);
                prho[facei] = volMixture_.rho(pp[facei], pT[facei]);

                palpha[facei] =
                    volMixture_.kappa(pp[facei], pT[facei])
                  / mixture_.Cpv(pp[facei], pT[facei]);
            }
        }
        else
        {
            forAll(pT, facei)
            {
                const typename MixtureType::thermoType& mixture_ =
                    this->patchFaceMixture(patchi, facei);

                const typename MixtureType::thermoType& volMixture_ =
                    this->patchFaceVolMixture
                    (
                        pp[facei],
                        pT[facei],
                        patchi,
                        facei
                    );

                pT[facei] = mixture_.THE(phe[facei], pp[facei] ,pT[facei]);
                prho[facei] = volMixture_.rho(pp[facei], pT[facei]);

                palpha[facei] =
                    volMixture_.kappa(pp[facei], pT[facei])
                  / mixture_.Cpv(pp[facei], pT[facei]);
            }
        }
    }
    
    this->alpha_.correctBoundaryConditions();
}


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

template<class BasicSolidThermo, class MixtureType>
Foam::heSolidThermo<BasicSolidThermo, MixtureType>::
heSolidThermo
(
    const fvMesh& mesh,
    const word& phaseName
)
:
    heThermo<BasicSolidThermo, MixtureType>(mesh, phaseName)
{
    calculate();
}


template<class BasicSolidThermo, class MixtureType>
Foam::heSolidThermo<BasicSolidThermo, MixtureType>::
heSolidThermo
(
    const fvMesh& mesh,
    const dictionary& dict,
    const word& phaseName
)
:
    heThermo<BasicSolidThermo, MixtureType>(mesh, dict, phaseName)
{
    calculate();
}


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

template<class BasicSolidThermo, class MixtureType>
Foam::heSolidThermo<BasicSolidThermo, MixtureType>::~heSolidThermo()
{}


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

template<class BasicSolidThermo, class MixtureType>
void Foam::heSolidThermo<BasicSolidThermo, MixtureType>::correct()
{
    if (debug)
    {
        InfoInFunction << endl;
    }

    calculate();

    if (debug)
    {
        Info<< "    Finished" << endl;
    }
}


template<class BasicSolidThermo, class MixtureType>
Foam::tmp<Foam::volVectorField>
Foam::heSolidThermo<BasicSolidThermo, MixtureType>::Kappa() const
{
    const fvMesh& mesh = this->T_.mesh();

    tmp<volVectorField> tKappa
    (
        new volVectorField
        (
            IOobject
            (
                "Kappa",
                mesh.time().timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            mesh,
            dimEnergy/dimTime/dimLength/dimTemperature
        )
    );

    volVectorField& Kappa = tKappa.ref();
    vectorField& KappaCells = Kappa.primitiveFieldRef();
    const scalarField& TCells = this->T_;
    const scalarField& pCells = this->p_;

    forAll(KappaCells, celli)
    {
        Kappa[celli] =
            this->cellVolMixture
            (
                pCells[celli],
                TCells[celli],
                celli
            ).Kappa(pCells[celli], TCells[celli]);
    }

    volVectorField::Boundary& KappaBf = Kappa.boundaryFieldRef();

    forAll(KappaBf, patchi)
    {
        vectorField& Kappap = KappaBf[patchi];
        const scalarField& pT = this->T_.boundaryField()[patchi];
        const scalarField& pp = this->p_.boundaryField()[patchi];

        forAll(Kappap, facei)
        {
            Kappap[facei] =
                this->patchFaceVolMixture
                (
                    pp[facei],
                    pT[facei],
                    patchi,
                    facei
                ).Kappa(pp[facei], pT[facei]);
        }
    }

    return tKappa;
}


template<class BasicSolidThermo, class MixtureType>
Foam::tmp<Foam::vectorField>
Foam::heSolidThermo<BasicSolidThermo, MixtureType>::Kappa
(
    const label patchi
) const
{
    const scalarField& pp = this->p_.boundaryField()[patchi];
    const scalarField& Tp = this->T_.boundaryField()[patchi];
    tmp<vectorField> tKappa(new vectorField(pp.size()));

    vectorField& Kappap = tKappa.ref();

    forAll(Tp, facei)
    {
        Kappap[facei] =
            this->patchFaceVolMixture
            (
                pp[facei],
                Tp[facei],
                patchi,
                facei
            ).Kappa(pp[facei], Tp[facei]);
    }

    return tKappa;
}


// ************************************************************************* //
heSolidThermo.C (8,236 bytes)   

henry

2017-10-23 10:02

manager   ~0008920

> Are the selection of the interpolation scheme for the diffusion coefficient and the specification of the diffusion coefficient at the boundaries really not related to each other?

Not necessarily, some kind of wall-function might be needed at the boundary to handle the particular coefficient distribution near the wall due to turbulence for example. It would be too restrictive to hard-code boundary condition treatment into the interpolation schemes.

I will test the patch and apply.

henry

2017-11-07 12:20

manager   ~0009004

In basicThermo.C why have you changed the constructor from uninitializing to initializing to 0?

- dimensionSet(1, -1, -1, 0, 0)
+ dimensionedScalar
+ (
+ "zero",
+ dimensionSet(1, -1, -1, 0, 0),
+ Zero
+ )

If initialization is essential it is unlikely that a value of 0 is appropriate and the constructor would need to be changed to MUST_READ.

bjoern.pfeiffelmann

2017-11-07 12:45

reporter   ~0009005

I thought it is more compatible to previous releases because the creation of an alpha dictionary is optional then

henry

2017-11-07 13:00

manager   ~0009006

Right, that makes sense.

henry

2017-11-07 13:56

manager   ~0009011

resolved by commit 204c6ee449ddf6eb416f648ccbb4c18d4ac04002

Issue History

Date Modified Username Field Change
2017-10-06 13:35 bjoern.pfeiffelmann New Issue
2017-10-06 13:35 bjoern.pfeiffelmann File Added: harmonicPatch.zip
2017-10-06 13:35 bjoern.pfeiffelmann Tag Attached: Heat transfer
2017-10-06 13:35 bjoern.pfeiffelmann Tag Attached: harmonic
2017-10-06 13:35 bjoern.pfeiffelmann Tag Attached: Surface Interpolation
2017-10-06 13:36 bjoern.pfeiffelmann File Added: harmonicPatch.jpg
2017-10-06 13:36 bjoern.pfeiffelmann File Added: OF_vs_ANALYTIC_T.png
2017-10-06 13:37 bjoern.pfeiffelmann File Added: OF_vs_ANALYTIC_wallHeatFlux.png
2017-10-06 15:04 henry Note Added: 0008825
2017-10-06 18:21 StephanG File Added: wall1D.tar.gz
2017-10-06 18:21 StephanG Note Added: 0008829
2017-10-06 18:31 henry Note Added: 0008830
2017-10-06 19:03 StephanG Note Added: 0008831
2017-10-06 20:12 henry Note Added: 0008832
2017-10-09 07:53 bjoern.pfeiffelmann Note Added: 0008833
2017-10-09 08:20 henry Note Added: 0008834
2017-10-09 08:51 henry Note Edited: 0008834
2017-10-09 09:48 bjoern.pfeiffelmann Note Added: 0008835
2017-10-09 10:04 henry Note Added: 0008836
2017-10-09 14:23 StephanG Note Added: 0008843
2017-10-09 14:28 henry Note Added: 0008844
2017-10-10 07:33 bjoern.pfeiffelmann Note Added: 0008846
2017-10-10 08:06 henry Note Added: 0008847
2017-10-10 13:46 bjoern.pfeiffelmann Note Added: 0008851
2017-10-13 09:31 henry Assigned To => henry
2017-10-13 09:31 henry Status new => closed
2017-10-13 09:31 henry Resolution open => suspended
2017-10-13 09:31 henry Note Added: 0008865
2017-10-20 12:27 bjoern.pfeiffelmann Status closed => feedback
2017-10-20 12:27 bjoern.pfeiffelmann Resolution suspended => reopened
2017-10-20 12:27 bjoern.pfeiffelmann Note Added: 0008907
2017-10-20 12:28 bjoern.pfeiffelmann File Added: wall1Dharmonic.tar.gz
2017-10-20 12:52 henry Note Added: 0008908
2017-10-20 15:23 bjoern.pfeiffelmann Note Added: 0008913
2017-10-20 15:23 bjoern.pfeiffelmann Status feedback => assigned
2017-10-20 15:27 henry Note Added: 0008914
2017-10-23 07:58 bjoern.pfeiffelmann Note Added: 0008919
2017-10-23 08:04 bjoern.pfeiffelmann File Added: basicThermo.C
2017-10-23 08:04 bjoern.pfeiffelmann File Added: heSolidThermo.C
2017-10-23 10:02 henry Note Added: 0008920
2017-11-07 12:20 henry Note Added: 0009004
2017-11-07 12:45 bjoern.pfeiffelmann Note Added: 0009005
2017-11-07 13:00 henry Note Added: 0009006
2017-11-07 13:56 henry Status assigned => resolved
2017-11-07 13:56 henry Resolution reopened => fixed
2017-11-07 13:56 henry Fixed in Version => dev
2017-11-07 13:56 henry Note Added: 0009011