[taban]

Dear foamers,
I really need your help to solve my problems in running openFoam 2.2.1. I am trying to solve a convection and radiation case, so far with little success. The case is a 2D room (4m*8m) with 2 inlet and 2 outlet for air flow. There are 2 heaters at uniform temperature (850 K) on the ceiling ,the geometry is attached (geometry.png)
For this, I use the bouyantBoussinesqPimpleFoam solver, fvDOM radiation model and RNGkEpsilon for turbulence model. This case is steady state. I am running unsteady solver for a longer physical time until it reaches steady state.
Some boundary conditions are defined as follows ...
T:

      internalField   uniform 263;
boundaryField
{
     heater
    {
        type            fixedValue;
        value           uniform 850;
    }
    inletL
    {
        type            fixedValue;
        value           uniform 263;
    }
   inletR
    {
        type            fixedValue;
        value           uniform 263;
    }
    outletL
    {
        type            inletOutlet;
        value           uniform 263;
        inletValue      uniform 263;
    }
outletR
    {
        type            inletOutlet;
        value           uniform 263;
        inletValue      uniform 263;
    }
    sidewalls
    {
     type            fixedValue;
     value           uniform 263;
    }
   
   ".*"
       {
        type            groovyBC;
        variables       "gradT=Qr/0.02334;";// air conductivity= 0.02334
        fractionExpression "0";
        gradientExpression "gradT";
        value           uniform 263;
    }
}
    }

U:

internalField   uniform (0 0 0);

boundaryField
{
   inletL
    {
        type            fixedValue;
        value           uniform (0.213 0 0);
    }
    inletR
    {
        type            fixedValue;
        value           uniform (-0.213 0 0);
    }
    outletL
    {
        type            pressureInletOutletVelocity;
        value           uniform (0 0 0);
        inletValue      uniform (0 0 0);
    }
outletR
    {
        type            pressureInletOutletVelocity;
        value           uniform (0 0 0);
        inletValue      uniform (0 0 0);
    }

    ".*"
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
}

IDefault

internalField   uniform 0;

boundaryField
{
    "heater"
    {
        type            greyDiffusiveRadiation;
        T               T;
        emissivityMode  lookup;
        emissivity      uniform 0.9;
        value           uniform 0;
    }

   ".*"
    {
        type            greyDiffusiveRadiation;
        T               T;
        emissivityMode  lookup;
        emissivity      uniform 0.5;
        value           uniform 0;
    }
}

controlDict

application     buoyantBoussinesqPimpleFoam;

startFrom       latestTime;

startTime       0;

stopAt          endTime;

endTime         50;

deltaT          0.001;

writeControl    timeStep;

writeInterval   1000;

purgeWrite      0;

writeFormat     ascii;

writePrecision  9;

writeCompression compressed;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;

adjustTimeStep  yes;

maxCo           0.1;

functions
{
   probes
   {
      type             probes;
      functionObjectLibs   ("libsampling.so");
      enabled              true;
      outputControl        outputTime; //timeStep;
      outputInterval       1;
      
      fields 
      (
         T
         U
      );

      probeLocations
      (
          (2 0.0 0.0)
          (4 0 0.0)
          (6 0 0.0)
      );
   }
}
 
libs ( "libOpenFOAM.so" "libgroovyBC.so" ) ;

fvSchemes:

ddtSchemes
{
    default         Euler;
}
gradSchemes
{
    default         Gauss linear;
}
divSchemes
{
    default         none;
    div(phi,U)      Gauss upwind;
    div(phi,T)      Gauss upwind;
    div(phi,k)      Gauss upwind;
    div(phi,epsilon) Gauss upwind;
    div(phi,R)      Gauss upwind;
    div(R)          Gauss linear;
    div((nuEff*dev(T(grad(U))))) Gauss linear;
    div(Ji,Ii_h)    Gauss linearUpwind grad(Ii_h);//ezafe shod
}
laplacianSchemes
{
    default         none;
    laplacian(nuEff,U) Gauss linear uncorrected;
    laplacian(Dp,p_rgh) Gauss linear uncorrected;
    laplacian(alphaEff,T) Gauss linear uncorrected;
    laplacian(DkEff,k) Gauss linear uncorrected;
    laplacian(DepsilonEff,epsilon) Gauss linear uncorrected;
    laplacian(DREff,R) Gauss linear uncorrected;
}
interpolationSchemes
{
    default         linear;
}
snGradSchemes
{
    default         uncorrected;
}
fluxRequired
{
    default         no;
    p_rgh;
}

fvSolution

solvers
{
    p_rgh
    {
        solver              GAMG;
        tolerance           1e-4;
        relTol              0;
        smoother            symGaussSeidel;
        cacheAgglomeration  true;
        nCellsInCoarsestLevel   10;
        agglomerator    faceAreaPair;
        mergeLevels         1;
        maxIter             1;
        nPreSweeps          0;
        nPostSweeps         1;
    }
    p_rghFinal
    {
        $p_rgh;
        relTol          0;
    }
    "(U|T|k|epsilon|R)"
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-4;//6
        relTol          0.1;
    }
    "(U|T|k|epsilon|R)Final"
    {
        $U;
        relTol          0;
    }
 Ii
    {
        solver              GAMG;
        tolerance           1e-4;
        relTol              0;
        smoother            symGaussSeidel;
        cacheAgglomeration  true;
        nCellsInCoarsestLevel   10;
        agglomerator    faceAreaPair;
        mergeLevels         1;
        maxIter             1;
        nPreSweeps          0;
        nPostSweeps         1;
    }
}
PIMPLE
{
    momentumPredictor no;
    nOuterCorrectors 1;
    nCorrectors     2;
    nNonOrthogonalCorrectors 0;
    pRefCell        0;
    pRefValue       0;
}
relaxationFactors
{
    fields
    {
    }
    equations
    {
        "(U|T|k|epsilon|R)" 1;
        "(U|T|k|epsilon|R)Final" 1;
    }
}

radiation properties:

radiation       on;

radiationModel  fvDOM;

fvDOMCoeffs
{
    nPhi    80;          // azimuthal angles in PI/2 on X-Y.(from Y to X)
    nTheta  2;          // polar angles in PI (from Z to X-Y plane)
    convergence 0.001;   // convergence criteria for radiation iteration
    maxIter 2;         // maximum number of iterations
    cacheDiv    true;   // cache the div of the RTE equation.
// NOTE: Caching div is "only" accurate if the upwind scheme is used in
// div(Ji,Ii_h)
}

// Number of flow iterations per radiation iteration
solverFreq 20;//5

absorptionEmissionModel constantAbsorptionEmission;

constantAbsorptionEmissionCoeffs
{
   absorptivity    absorptivity    [ 0 -1 0 0 0 0 0 ] 0;
   emissivity      emissivity      [ 0 -1 0 0 0 0 0 ] 0;
   E               E   [ 1 -1 -3 0 0 0 0 ] 0;
}

scatterModel    none;

The correct solution isn't available but temperature contour and speed contour obviously shows that solution is incorrect (both images are attached). It would be nice if anyone can help me
1) I want the floor and cieling of the room define as adiabatic boundary so i use groovy bc as you see above in T code. Is it correct?
2) if i choose maxcorunt number > 0.1 in controlDict the solver crash after some timesteps. what is the problem?
3) Maybe i choose an improper scheme or other parameters in controldict , fvSolution or ..., please check them and tell me What might be the mistake? If it is required more information in details feel free to ask.
thanks in advance
Best regards