View Issue Details

IDProjectCategoryView StatusLast Update
0001704OpenFOAMBugpublic2015-05-25 17:36
Reporterzfaraday Assigned Tohenry  
PriorityhighSeveritymajorReproducibilityalways
Status resolvedResolutionfixed 
PlatformGNU/LinuxOSOpenSuSEOS Version12.3
Summary0001704: externalWallHeatFluxTemperature BC makes the solver crash when Qr is added to the balance
DescriptionAfter trying to solve some cases with buoyantSimpleFoam and chtMultiRegionSImpleFOam using this BC with the option of adding Qr to the balance I never got any other solution than a crash. After an intensive analysis (http://www.cfd-online.com/Forums/openfoam-solving/143632-problems-adding-qr-field-externalwallheatfluxtemperature-bc.html#post547223) I found out that the problem was that when Qr is added to the balance it makes refValue() get negative values in the first iterations.

In order to solve this problem I tried to implement a relaxation factor for Qr like it is done in thermalBaffle1D BC with hopeful results. The solver converged with a relaxation factor of 0.1.

I attached a case where the problem can be tested, it is a modified version of the tutorial case "circuitBoardCooling". Here (http://www.cfd-online.com/Forums/openfoam-bugs/153082-thermalbaffle1d-showing-bad-behavior.html#post546783) its changes are explained in order to see the crash.
Steps To Reproduceexecute ./Allrun script

The solver should hopefully crash after 16th timestep.
TagsNo tags attached.

Activities

zfaraday

2015-05-22 13:24

reporter  

zfaraday

2015-05-22 13:25

reporter  

myExternalWallHeatFluxTemperatureFvPatchScalarField.C (11,979 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2014 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 "myExternalWallHeatFluxTemperatureFvPatchScalarField.H"
#include "addToRunTimeSelectionTable.H"
#include "fvPatchFieldMapper.H"
#include "volFields.H"
#include "mappedPatchBase.H"

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

namespace Foam
{

template<>
const char*
NamedEnum
<myExternalWallHeatFluxTemperatureFvPatchScalarField::operationMode, 3>::names[]=
{
    "fixed_heat_flux",
    "fixed_heat_transfer_coefficient",
    "unknown"
};

const NamedEnum
<
    myExternalWallHeatFluxTemperatureFvPatchScalarField::operationMode, 3
>
myExternalWallHeatFluxTemperatureFvPatchScalarField::operationModeNames;

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

} // End namespace Foam

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

Foam::myExternalWallHeatFluxTemperatureFvPatchScalarField::
myExternalWallHeatFluxTemperatureFvPatchScalarField
(
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF
)
:
    mixedFvPatchScalarField(p, iF),
    temperatureCoupledBase(patch(), "undefined", "undefined", "undefined-K"),
    mode_(unknown),
    q_(p.size(), 0.0),
    h_(p.size(), 0.0),
    Ta_(p.size(), 0.0),
    QrName_("undefined-Qr"),
    thicknessLayers_(),
    kappaLayers_(),
    //Fields added
    QrRelaxation_(1),
    QrPrevious_(p.size())
{
    refValue() = 0.0;
    refGrad() = 0.0;
    valueFraction() = 1.0;
}


Foam::myExternalWallHeatFluxTemperatureFvPatchScalarField::
myExternalWallHeatFluxTemperatureFvPatchScalarField
(
    const myExternalWallHeatFluxTemperatureFvPatchScalarField& ptf,
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF,
    const fvPatchFieldMapper& mapper
)
:
    mixedFvPatchScalarField(ptf, p, iF, mapper),
    temperatureCoupledBase(patch(), ptf),
    mode_(ptf.mode_),
    q_(ptf.q_, mapper),
    h_(ptf.h_, mapper),
    Ta_(ptf.Ta_, mapper),
    QrName_(ptf.QrName_),
    thicknessLayers_(ptf.thicknessLayers_),
    kappaLayers_(ptf.kappaLayers_),
    //fields added
    QrRelaxation_(ptf.QrRelaxation_),
    QrPrevious_(ptf.QrPrevious_, mapper)
{}


Foam::myExternalWallHeatFluxTemperatureFvPatchScalarField::
myExternalWallHeatFluxTemperatureFvPatchScalarField
(
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF,
    const dictionary& dict
)
:
    mixedFvPatchScalarField(p, iF),
    temperatureCoupledBase(patch(), dict),
    mode_(unknown),
    q_(p.size(), 0.0),
    h_(p.size(), 0.0),
    Ta_(p.size(), 0.0),
    QrName_(dict.lookupOrDefault<word>("Qr", "none")),
    thicknessLayers_(),
    kappaLayers_(),
    //fields added
    QrRelaxation_(dict.lookupOrDefault<scalar>("relaxation", 1)),
    QrPrevious_(p.size(), 0.0)
{
    if (dict.found("q") && !dict.found("h") && !dict.found("Ta"))
    {
        mode_ = fixedHeatFlux;
        q_ = scalarField("q", dict, p.size());
    }
    else if (dict.found("h") && dict.found("Ta") && !dict.found("q"))
    {
        mode_ = fixedHeatTransferCoeff;
        h_ = scalarField("h", dict, p.size());
        Ta_ = scalarField("Ta", dict, p.size());
        if (dict.found("thicknessLayers"))
        {
            dict.lookup("thicknessLayers") >> thicknessLayers_;
            dict.lookup("kappaLayers") >> kappaLayers_;
        }
    }
    else
    {
        FatalErrorIn
        (
            "myExternalWallHeatFluxTemperatureFvPatchScalarField::"
            "myExternalWallHeatFluxTemperatureFvPatchScalarField\n"
            "(\n"
            "    const fvPatch& p,\n"
            "    const DimensionedField<scalar, volMesh>& iF,\n"
            "    const dictionary& dict\n"
            ")\n"
        )   << "\n patch type '" << p.type()
            << "' either q or h and Ta were not found '"
            << "\n for patch " << p.name()
            << " of field " << dimensionedInternalField().name()
            << " in file " << dimensionedInternalField().objectPath()
            << exit(FatalError);
    }

    fvPatchScalarField::operator=(scalarField("value", dict, p.size()));
    //Added
    if (dict.found("QrPrevious"))
    {
        QrPrevious_ = scalarField("QrPrevious", dict, p.size());
    }
//
    if (dict.found("refValue"))
    {
        // Full restart
        refValue() = scalarField("refValue", dict, p.size());
        refGrad() = scalarField("refGradient", dict, p.size());
        valueFraction() = scalarField("valueFraction", dict, p.size());
    }
    else
    {
        // Start from user entered data. Assume fixedValue.
        refValue() = *this;
        refGrad() = 0.0;
        valueFraction() = 1.0;
    }
}


Foam::myExternalWallHeatFluxTemperatureFvPatchScalarField::
myExternalWallHeatFluxTemperatureFvPatchScalarField
(
    const myExternalWallHeatFluxTemperatureFvPatchScalarField& tppsf
)
:
    mixedFvPatchScalarField(tppsf),
    temperatureCoupledBase(tppsf),
    mode_(tppsf.mode_),
    q_(tppsf.q_),
    h_(tppsf.h_),
    Ta_(tppsf.Ta_),
    QrName_(tppsf.QrName_),
    thicknessLayers_(tppsf.thicknessLayers_),
    kappaLayers_(tppsf.kappaLayers_),
    //field added
    QrRelaxation_(tppsf.QrRelaxation_),
    QrPrevious_(tppsf.QrPrevious_)
{}


Foam::myExternalWallHeatFluxTemperatureFvPatchScalarField::
myExternalWallHeatFluxTemperatureFvPatchScalarField
(
    const myExternalWallHeatFluxTemperatureFvPatchScalarField& tppsf,
    const DimensionedField<scalar, volMesh>& iF
)
:
    mixedFvPatchScalarField(tppsf, iF),
    temperatureCoupledBase(patch(), tppsf),
    mode_(tppsf.mode_),
    q_(tppsf.q_),
    h_(tppsf.h_),
    Ta_(tppsf.Ta_),
    QrName_(tppsf.QrName_),
    thicknessLayers_(tppsf.thicknessLayers_),
    kappaLayers_(tppsf.kappaLayers_),
    //fields added
    QrRelaxation_(tppsf.QrRelaxation_),
    QrPrevious_(tppsf.QrPrevious_)
{}


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

void Foam::myExternalWallHeatFluxTemperatureFvPatchScalarField::autoMap
(
    const fvPatchFieldMapper& m
)
{
    mixedFvPatchScalarField::autoMap(m);
    q_.autoMap(m);
    h_.autoMap(m);
    Ta_.autoMap(m);
}


void Foam::myExternalWallHeatFluxTemperatureFvPatchScalarField::rmap
(
    const fvPatchScalarField& ptf,
    const labelList& addr
)
{
    mixedFvPatchScalarField::rmap(ptf, addr);

    const myExternalWallHeatFluxTemperatureFvPatchScalarField& tiptf =
        refCast<const myExternalWallHeatFluxTemperatureFvPatchScalarField>(ptf);

    q_.rmap(tiptf.q_, addr);
    h_.rmap(tiptf.h_, addr);
    Ta_.rmap(tiptf.Ta_, addr);
}


void Foam::myExternalWallHeatFluxTemperatureFvPatchScalarField::updateCoeffs()
{
    if (updated())
    {
        return;
    }

    const scalarField Tp(*this);
    scalarField hp(patch().size(), 0.0);

    scalarField Qr(Tp.size(), 0.0);
    if (QrName_ != "none")
    {
        Qr = patch().lookupPatchField<volScalarField, scalar>(QrName_);
        //Relaxation of Qr like in thermalBaffle1D
        Qr = QrRelaxation_*Qr + (1.0 - QrRelaxation_)*QrPrevious_;
        QrPrevious_ = Qr;
    }

    switch (mode_)
    {
        case fixedHeatFlux:
        {
            break;
        }
        case fixedHeatTransferCoeff:
        {
            scalar totalSolidRes = 0.0;
            if (thicknessLayers_.size() > 0)
            {
                forAll (thicknessLayers_, iLayer)
                {
                    const scalar l = thicknessLayers_[iLayer];
                    if (kappaLayers_[iLayer] > 0.0)
                    {
                        totalSolidRes += l/kappaLayers_[iLayer];
                    }
                }
            }
            hp = 1.0/(1.0/h_ + totalSolidRes);
            break;
        }
        default:
        {
            FatalErrorIn
            (
                "myExternalWallHeatFluxTemperatureFvPatchScalarField"
                "::updateCoeffs()"
            )   << "Illegal heat flux mode " << operationModeNames[mode_]
                << exit(FatalError);
        }
    }

    if (mode_ == fixedHeatFlux)
    {
        refGrad() =  (q_ + Qr)/kappa(Tp);
        refValue() =  0.0;
        valueFraction() = 0.0;
    }
    else if (mode_ == fixedHeatTransferCoeff)
    {
        //control outputs to check possible errors
/*        Info << "patchName: " << patch().name() << endl;
        Info << "patchType: " << patch().type() << endl;
        Info << "Qr: " << Qr << endl;
*/
//        Qr /= Tp;
        refGrad() =  0.0;
        refValue() =  Ta_ + Qr/hp;
        valueFraction() = hp / (hp + kappa(Tp)*patch().deltaCoeffs());

        //control outputs to check possible errors
/*        Info << "hp: " << hp << endl;
        Info << "Qr/Tp: " << Qr << endl;
        Info << "Tp: " << Tp << endl;
        Info << "KDelta: " << kappa(Tp)*patch().deltaCoeffs() << endl;
        Info << "refValue(): " << refValue() << endl;
        Info << "valueFraction(): " << valueFraction() << endl;
*/
    }

    mixedFvPatchScalarField::updateCoeffs();

    if (debug)
    {
        scalar Q = gSum(kappa(Tp)*patch().magSf()*snGrad());

        Info<< patch().boundaryMesh().mesh().name() << ':'
            << patch().name() << ':'
            << this->dimensionedInternalField().name() << " :"
            << " heat transfer rate:" << Q
            << " walltemperature "
            << " min:" << gMin(*this)
            << " max:" << gMax(*this)
            << " avg:" << gAverage(*this)
            << endl;
    }
}


void Foam::myExternalWallHeatFluxTemperatureFvPatchScalarField::write
(
    Ostream& os
) const
{
    mixedFvPatchScalarField::write(os);
    temperatureCoupledBase::write(os);

    QrPrevious_.writeEntry("QrPrevious", os);//
    os.writeKeyword("Qr")<< QrName_ << token::END_STATEMENT << nl;
    os.writeKeyword("relaxation")<< QrRelaxation_
        << token::END_STATEMENT << nl;//

    switch (mode_)
    {

        case fixedHeatFlux:
        {
            q_.writeEntry("q", os);
            break;
        }
        case fixedHeatTransferCoeff:
        {
            h_.writeEntry("h", os);
            Ta_.writeEntry("Ta", os);
            thicknessLayers_.writeEntry("thicknessLayers", os);
            kappaLayers_.writeEntry("kappaLayers", os);
            break;
        }
        default:
        {
            FatalErrorIn
            (
                "void myExternalWallHeatFluxTemperatureFvPatchScalarField::write"
                "("
                    "Ostream& os"
                ") const"
            )   << "Illegal heat flux mode " << operationModeNames[mode_]
                << abort(FatalError);
        }
    }
}


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

namespace Foam
{
    makePatchTypeField
    (
        fvPatchScalarField,
        myExternalWallHeatFluxTemperatureFvPatchScalarField
    );
}

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

zfaraday

2015-05-22 13:31

reporter   ~0004785

I forgot to mention that Qr option is not activated by default in the case attached. It must be activated in the "0.org/T" file.

henry

2015-05-24 19:39

manager   ~0004792

I have added relaxation of the radiative contribution to the externalWallHeatFluxTemperature BC as requested but I am unable to run your case with or without this change:

Starting time loop

Time = 1

DILUPBiCG: Solving for h, Initial residual = nan, Final residual = nan, No Iterations 1001
Radiation solver iter: 0

Also your case has a lot of external dependencies which are not provided with OpenFOAM which makes it difficult to test.

Please test the change yourself to ensure the relaxation operates as you expect.

commit 24059d8139ef946984e7f262ba1272366e4daf5f in OpenFOAM-2.4.x

zfaraday

2015-05-24 21:31

reporter   ~0004794

Sorry Henry, but I don't understand what you mean with external dependencies, I only took a tutorial case and modified its boundary conditions using other ones provided with OpenFOAM. Besides that I just realized that I used some functions and libraries in controlDict that maybe are causing your troubles. This is the only thing I think that can cause a trouble... I apologize for that, I only added these functions to check the energy balance and I forgot to remove it when I uploaded here.

I see you only applied the change in OF v2.4.x , may you apply it to v2.3.x too? This the version I'm using at the moment and I'm not going to install 2.4.x for now.

Btw, I just checked the modifications and I think that some lines are wrong, these are lines 219-220 and this is how I think they should look:

     q_(tppsf.q_),
     h_(tppsf.h_),
     Ta_(tppsf.Ta_),
+ QrPrevious_(ptf.QrPrevious_), ---> QrPrevious_(tppsf.QrPrevious_), //219
+ QrRelaxation_(ptf.QrRelaxation_), ---> QrRelaxation_(tppsf.QrRelaxation_), //220
     QrName_(tppsf.QrName_),
     thicknessLayers_(tppsf.thicknessLayers_),
     kappaLayers_(tppsf.kappaLayers_)

Many thanks for applying the changes in order to solve this issue.

Best regards,

Alex

wyldckat

2015-05-24 21:36

updater   ~0004795

Just a quick note on how to apply the relevant commit on 2.3.x manually:

  foam

  wget https://github.com/OpenFOAM/OpenFOAM-2.4.x/commit/24059d8139ef946984e7f262ba1272366e4daf5f.patch

  patch -p1 < 24059d8139ef946984e7f262ba1272366e4daf5f.patch

  wmake src/turbulenceModels/compressible/turbulenceModel

henry

2015-05-24 22:21

manager   ~0004796

Sorry about the compilation error; juggling too many things. I have pushed the fix to OpenFOAM-2.4.x. I will make the same change to 2.3.x if it works as expected.

zfaraday

2015-05-25 17:28

reporter   ~0004808

Thanks for the tip Bruno!

Now I just tested the new patch specification and it works exactly as expected in the test case attached.

Bug resolved! Thanks Henry

henry

2015-05-25 17:36

manager   ~0004809

Resolved by commit 95a09ca08c5d5e1b9797893df72451d8a8d41317

Issue History

Date Modified Username Field Change
2015-05-22 13:24 zfaraday New Issue
2015-05-22 13:24 zfaraday File Added: circuitBoardCooling_rad.zip
2015-05-22 13:25 zfaraday File Added: myExternalWallHeatFluxTemperatureFvPatchScalarField.C
2015-05-22 13:31 zfaraday Note Added: 0004785
2015-05-24 19:39 henry Note Added: 0004792
2015-05-24 21:31 zfaraday Note Added: 0004794
2015-05-24 21:36 wyldckat Note Added: 0004795
2015-05-24 22:21 henry Note Added: 0004796
2015-05-25 17:28 zfaraday Note Added: 0004808
2015-05-25 17:36 henry Note Added: 0004809
2015-05-25 17:36 henry Status new => resolved
2015-05-25 17:36 henry Resolution open => fixed
2015-05-25 17:36 henry Assigned To => henry