View Issue Details

Jump to NotesJump to History
IDProjectCategoryView StatusDate SubmittedLast Update
0001369OpenFOAMBugpublic2014-08-07 14:162014-12-29 11:09
Reporteruser700Assigned Tohenry  
PrioritynormalSeveritymajorReproducibilityalways
Status resolvedResolutionfixed 
Platformany 
Summary0001369: turbulentTemperatureRadCoupledMixedFvPatchScalarField::updateCoeffs(): refValue can become negative
DescriptionIn turbulentTemperatureRadCoupledMixedFvPatchScalarField::updateCoeffs()
a field alpha is computed as KDeltaNbr - (Qr + QrNbr)/Tp, which can become negative, if Qr is large, and results in a negative refValue later. This is passed to a thermo.He() function as a temperature in enthalpy based solver, leading to a crash. (Or to nonsensical calculations of enthalpy as function of temperature).
I propose to include the radiative contribution in the refGrad instead, as refGrad=(Qr+QrNbr)/kappa. I did this and my case where I ran into this problem is working fine now.
TagsNo tags attached.

Activities

user700

2014-08-08 12:33

  ~0003199

I forgot to mention: The trouble really arises with temperature dependent heat capacities ...

user21

2014-08-11 09:23

  ~0003201

Could you provide a simple case for testing?
Thanks
Sergio

user700

2014-08-11 11:06

  ~0003202

Might take a week because I wont be here but I will try to set up a case based on the tutorial files.

user700

2014-08-13 14:29

  ~0003203

Found an easy way to reproduce this:

I am using OpenFoam 2.2.x, but as far as I see it should be the same with 2.3.

1.) Take tutorial case multiRegionHeaterRadiation
2.) change constant/heater/thermophysicalProperties to this:

thermoType
{
    type heSolidThermo;
    mixture pureMixture;
    //transport constIso;
    transport polynomial;
    //thermo hConst;
    thermo hPolynomial;
    equationOfState rhoConst;
    specie specie;
    energy sensibleEnthalpy;
}

mixture
{
    specie
    {
        nMoles 1;
        molWeight 12;
    }

    transport
    {
        //kappa 80;
    kappaCoeffs<8> (80 0 0 0 0 0 0 0);
    }

    thermodynamics
    {
        Hf 0;
    Sf 0.;
        // Cp 450;
    CpCoeffs<8> (200 0.6 0 0 0 0 0 0);
    }

    equationOfState
    {
        rho 8000;
    }
}

3.) change the boundary condition in 0/heater/T for minY from fixedValue 300 to fixedValue 1500

4.) run chtMultiRegionSimpleFoam and at t=374 this will (hopefully :) happen:
--> FOAM FATAL ERROR:
Maximum number of iterations exceeded

    From function thermo<Thermo, Type>::T(scalar f, scalar T0, scalar (thermo<Thermo, Type>::*F)(const scalar) const, scalar (thermo<Thermo, Type>::*dFdT)(const scalar) const, scalar (thermo<Thermo, Type>::*limit)(const scalar) const) const
  
5.) change in updateCoeffs() to:
    valueFraction() = KDeltaNbr/(KDeltaNbr+KDelta);
    refValue() = TcNbr;
    refGrad() = (Qr+QrNbr)/kappa(*this);

run again, for me it works, then.

user700

2014-08-26 12:14

  ~0003220

To be a little more precise: I think here Tw.refValue() should really evaluate to a temperature, not some effective value, because it is passed to the thermo to calculate the enthalpy:

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

    const basicThermo& thermo = basicThermo::lookupThermo(*this);
    const label patchi = patch().index();

    const scalarField& pw = thermo.p().boundaryField()[patchi];
    mixedFvPatchScalarField& Tw = refCast<mixedFvPatchScalarField>
    (
        const_cast<fvPatchScalarField&>(thermo.T().boundaryField()[patchi])
    );

    Tw.evaluate();

    valueFraction() = Tw.valueFraction();
    refValue() = thermo.he(pw, Tw.refValue(), patchi);
    refGrad() =
        thermo.Cpv(pw, Tw, patchi)*Tw.refGrad()
      + patch().deltaCoeffs()*
        (
            thermo.he(pw, Tw, patchi)
          - thermo.he(pw, Tw, patch().faceCells())
        );

    mixedFvPatchScalarField::updateCoeffs();
}

henry

2014-12-29 11:09

manager   ~0003375

Thanks for the bug-report and analysis of the problem. I have now studied the issue carefully and reproduced the problem in OpenFOAM-2.3.x based on your suggestion but I had to significantly increase the heater temperature to cause a crash.

As you suggest there are two actually two issues, one relating to a potential failure if Qr+QrNbr causes the reference temperature to become negative and the second relating to the local-linearization used in the T->h->T boundary condition conversions which requires the reference temperature to be close to the boundary temperature in the case of non-linear properties. Both of these issues are resolved by your suggested fix:

    valueFraction() = KDeltaNbr/(KDeltaNbr+KDelta);
    refValue() = TcNbr;
    refGrad() = (Qr+QrNbr)/kappa(*this);

which I have applied to OpenFOAM-2.3.x and OpenFOAM-dev.

Resolved by commit a01b5bb4a064656d07825359784fa34424c2b894

Issue History

Date Modified Username Field Change
2014-08-07 14:16 user700 New Issue
2014-08-08 12:33 user700 Note Added: 0003199
2014-08-11 09:23 user21 Note Added: 0003201
2014-08-11 11:06 user700 Note Added: 0003202
2014-08-13 14:29 user700 Note Added: 0003203
2014-08-26 12:14 user700 Note Added: 0003220
2014-12-29 11:09 henry Note Added: 0003375
2014-12-29 11:09 henry Status new => resolved
2014-12-29 11:09 henry Resolution open => fixed
2014-12-29 11:09 henry Assigned To => henry