View Issue Details

IDProjectCategoryView StatusLast Update
0000806OpenFOAMBugpublic2013-04-05 16:31
Reporterjherb Assigned Tohenry  
PrioritynormalSeveritycrashReproducibilityalways
Status closedResolutionfixed 
PlatformLinuxOSOpenSuseOS Version11.3
Summary0000806: turbulentHeatFluxTemperature causes crash if yPlus is too low
DescriptionUsing the solver buoyantBoussinesqSimpleFoam and the boundary condition turbulentHeatFluxTemperature at a wall causes a crash in the solver for T if the yPlus value is too low at this wall.

The reason is, that the gradient of T at that wall is calculated in OpenFOAM-2.2.x/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/turbulentHeatFluxTemperature/turbulentHeatFluxTemperatureFvPatchScalarField.C:
gradient() = q_/(Ap*Cp0*alphaEffp);
or
gradient() = q_/(Cp0*alphaEffp);

But alphaEffp is set in /OpenFOAM-2.2.x/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H as:
kappat = turbulence->nut()/Prt;
kappat.correctBoundaryConditions();

nut() is calculated at the wall by the wall function, e.g. in /OpenFOAM-2.2.x/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkWallFunction/nutkWallFunctionFvPatchScalarField.C:
        if (yPlus > yPlusLam_)
        {
            nutw[faceI] = nuw[faceI]*(yPlus*kappa_/log(E_*yPlus) - 1.0);
        }

Now if yPlus is too low, nut is 0 so the equation for the gradient divides by 0.

And obvious solution would be to use a coarser grid, but this does not make sense (the SST turbulence model is used, so there should be no lower limit on yPlus).

How to fix this problem?
TagsNo tags attached.

Activities

henry

2013-04-04 20:01

manager   ~0002088

In TEqn.H

volScalarField kappaEff("kappaEff", turbulence->nu()/Pr + kappat);

are you saying that turbulence->nu() is also 0?

jherb

2013-04-04 21:41

reporter   ~0002089

I patched turbulentHeatFluxTemperatureFvPatchScalarField.C:
        case hsFlux:
        {
            Info << alphaEffp << endl;
            gradient() = q_/(Cp0*alphaEffp);
            break;
        }

This is the output (with three heated walls):
Starting time loop

Time = 1

DILUPBiCG: Solving for Ux, Initial residual = 1, Final residual = 0.0103772, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 1, Final residual = 0.0103764, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 1, Final residual = 0.00321814, No Iterations 1
587{0}
125{0}
82{0}
DILUPBiCG: Solving for T, Initial residual = nan, Final residual = nan, No Iterations 1001
DICPCG: Solving for p_rgh, Initial residual = nan, Final residual = nan, No Iterations 1001
time step continuity errors : sum local = nan, global = -nan, cumulative = -nan

So yes, it looks like nu is also 0.

jherb

2013-04-04 21:45

reporter   ~0002090

Ok, now I'm confused:
If I add this:
    const scalarField& nuw =
        patch().lookupPatchField<volScalarField, scalar>("nu");

    // retrieve (constant) specific heat capacity from transport dictionary
    const turbulenceModel& turbulence =
        db().lookupObject<turbulenceModel>("turbulenceModel");
    const scalar Cp0(readScalar(turbulence.transport().lookup("Cp0")));

    switch (heatSource_)
    {
        case hsPower:
        {
            const scalar Ap = gSum(patch().magSf());
            gradient() = q_/(Ap*Cp0*alphaEffp);
            break;
        }
        case hsFlux:
        {
            Info << alphaEffp << endl;
            Info << nuw << endl;
            gradient() = q_/(Cp0*alphaEffp);
            break;
        }

The output is:
Time = 1

DILUPBiCG: Solving for Ux, Initial residual = 1, Final residual = 0.0103772, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 1, Final residual = 0.0103764, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 1, Final residual = 0.00321814, No Iterations 1
587{0}
587{6.53e-05}
125{0}
125{6.53e-05}
82{0}
82{6.53e-05}
DILUPBiCG: Solving for T, Initial residual = nan, Final residual = nan, No Iterations 1001

So nu is not 0 at the wall but the line
    volScalarField kappaEff("kappaEff", turbulence->nu()/Pr + kappat);
seems not to have any effect?!?

jherb

2013-04-05 10:28

reporter   ~0002091

I used the wrong boundary condition setup in 0/T:
        alphaEff kappat;
instead of
        alphaEff kappaEff;

As there is no tutorial for this boundary condition I probably took my setup from the internet (perhaps: http://www.cfd-online.com/Forums/openfoam-solving/114418-buoyantboussinesqpimplefoam-turbulentheatfluxtemperature.html)

So this bug report can be closed.

Issue History

Date Modified Username Field Change
2013-04-04 19:21 jherb New Issue
2013-04-04 20:01 henry Note Added: 0002088
2013-04-04 21:41 jherb Note Added: 0002089
2013-04-04 21:45 jherb Note Added: 0002090
2013-04-05 10:28 jherb Note Added: 0002091
2013-04-05 16:31 henry Status new => closed
2013-04-05 16:31 henry Assigned To => henry
2013-04-05 16:31 henry Resolution open => fixed