View Issue Details Jump to Notes ] Issue History ] Print ]
IDProjectCategoryView StatusDate SubmittedLast Update
0000806OpenFOAM[All Projects] Bugpublic2013-04-04 18:212013-04-05 15:31
Reporterjherb 
Assigned Tohenry 
PrioritynormalSeveritycrashReproducibilityalways
StatusclosedResolutionfixed 
PlatformLinuxOSOpenSuseOS Version11.3
Product Version2.2.x 
Target VersionFixed in Version 
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.
Attached Files

- Relationships

-  Notes
(0002088)
henry (manager)
2013-04-04 19:01

In TEqn.H

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

are you saying that turbulence->nu() is also 0?
(0002089)
jherb (reporter)
2013-04-04 20:41

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.
(0002090)
jherb (reporter)
2013-04-04 20:45

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?!?
(0002091)
jherb (reporter)
2013-04-05 09:28

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 18:21 jherb New Issue
2013-04-04 19:01 henry Note Added: 0002088
2013-04-04 20:41 jherb Note Added: 0002089
2013-04-04 20:45 jherb Note Added: 0002090
2013-04-05 09:28 jherb Note Added: 0002091
2013-04-05 15:31 henry Status new => closed
2013-04-05 15:31 henry Assigned To => henry
2013-04-05 15:31 henry Resolution open => fixed