CFD Online Discussion Forums - Problems adding Qr field in externalWallHeatFluxTemperature BC

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)

- - Problems adding Qr field in externalWallHeatFluxTemperature BC (https://www.cfd-online.com/Forums/openfoam-solving/143632-problems-adding-qr-field-externalwallheatfluxtemperature-bc.html)

Problems adding Qr field in externalWallHeatFluxTemperature BC

Dear foamers,

Recently I found out a new feature added to OF 2.3.x for chtMultiRegion* solvers. This is the frozenFlow option that allows to only develop the energy equation in a fluid region. That, I think, permits to save much time when solving cases in vacuum conditions where pressure and velocity equations are not needed at all. I have been working for quite long in cases in vacuum conditions where only radiative heat transfer is important, so when I discovered this feature I quickly updated OF in order to try it. While I was waiting for the compilation to finish I kept looking for new features I still didn't know to be used in my old cases of heat transfer (radiative mainly). I got glad when I discovered that now the externalWallHeatFluxTemperature boundary condition allows to include Qr in the computation of T at the boundaries, it's something I used to miss some months ago...

Once OF finished its compilation I got back to an old simple test case I used to work in. It consist of a solid object and a heater (created as a thermal baffle) that heats up all the domain. Both are surrounded by a fluid domain (vacuum in this case). The geometry can be seen in the attached figure.

All external boundaries of the air domain are defined as externalWallHeatFluxTemperature. All of them, except one, have a fixed heat flux (q=0, isolated). The only one that isn't isolated (shown in grey in the picture) is defined as follows:

Code:

     maxY 

     { 

          type            externalWallHeatFluxTemperature; 

          kappa           fluidThermo;       // fluidThermo, solidThermo or 

                                         // lookup 

          Ta              uniform 300.0;     // ambient temperature /[K] 

          h               uniform 1.0;      // heat transfer coeff /[W/K/m2] 

          thicknessLayers (0.02); 

          kappaLayers     (35); 

          value           uniform 600;     // initial temperature / [K] 

          kappaName       none; 

          Qr              Qr; 

      }

I also added the Qr field to all other external walls.

So far it's the same procedure I used before, the only difference is the addition of the Qr field at the boundaries and the use of the frozenFlow option. However, when I run the solver it blows up at the 3rd iteration everytime. I attach my fvSolution and fvSchemes files for the air region so that you can evaluate them and see if I am not using the proper implementation.

fvSchemes

Code:

/*--------------------------------*- C++ -*----------------------------------*\

| =========                 |                                                 |

| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |

|  \\    /   O peration     | Version:  2.2.1                                 |

|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |

|    \\/     M anipulation  |                                                 |

\*---------------------------------------------------------------------------*/

FoamFile

{

    version     2.0;

    format      ascii;

    class       dictionary;

    object      fvSchemes;

}

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



ddtSchemes

{

    default steadyState;

}



gradSchemes

{

    default         Gauss linear;//cellLimited Gauss linear 0.5; //Gauss linear;

}  



divSchemes

{

    default         none;

    div(phi,U)      bounded Gauss upwind;

    div(phi,K)      bounded Gauss upwind;

    div(phi,h)      bounded Gauss upwind;

    div(phi,k)      bounded Gauss upwind;

    div(phi,K)      bounded Gauss upwind;

    div(phi,epsilon) bounded Gauss upwind;

    div(phi,R)      bounded Gauss upwind;

    div(R)          Gauss linear;

    div(Ji,Ii_h)    bounded Gauss linearUpwind grad(U);

    div((muEff*dev2(T(grad(U))))) Gauss linear;

}



laplacianSchemes

{

    default         none;

    laplacian(muEff,U) Gauss linear uncorrected;

    laplacian(Dp,p_rgh) Gauss linear uncorrected;

    laplacian(alphaEff,h) Gauss linear uncorrected;

    laplacian(DkEff,k) Gauss linear uncorrected;

    laplacian(DepsilonEff,epsilon) Gauss linear uncorrected;

    laplacian(DREff,R) Gauss linear uncorrected;

    laplacian(gammaRad,G) Gauss linear uncorrected;

    laplacian(rhorAUf,p_rgh) Gauss linear uncorrected;

}



interpolationSchemes

{

    default         linear;

}



snGradSchemes

{

    default         uncorrected;

}



fluxRequired

{

    default         no;

    p_rgh;

}



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

fvSolution

Code:

/*--------------------------------*- C++ -*----------------------------------*\

| =========                 |                                                 |

| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |

|  \\    /   O peration     | Version:  2.2.1                                 |

|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |

|    \\/     M anipulation  |                                                 |

\*---------------------------------------------------------------------------*/

FoamFile

{

    version     2.0;

    format      ascii;

    class       dictionary;

    object      fvSolution;

}

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



solvers

{

    h

    {

        solver           PBiCG;

        preconditioner   DILU;

        tolerance        1e-5;//1e-7;

        relTol           0;//0.01;

    }



    Ii

    {

        solver              GAMG;

        tolerance           1e-4;

        relTol              0;

        smoother            symGaussSeidel;

        cacheAgglomeration  true;

        nCellsInCoarsestLevel   10;

        agglomerator        faceAreaPair;

        mergeLevels         1;

        maxIter             5;

        nPreSweeps          0;

        nPostSweeps         1;

    }

/*    {

        $h;

        tolerance        1e-4;//1e-7;

        relTol           0.01;

    }

*/

}



SIMPLE

{

    momentumPredictor off;

    nNonOrthogonalCorrectors 0;

    pRefCell        0;

    pRefValue       1e-2;//100000;

    rhoMin          rhoMin [1 -3 0 0 0] 1e-10;

    rhoMax          rhoMax [1 -3 0 0 0] 20;

    frozenFlow      true;

    residualControl

    {

/*        h               1e-3;

        G               1e-3;

        "ILambda.*"     1e-3;

*/    }

}



relaxationFactors

{

    fields

    {

    }

    equations

    {

        h               0.3;//0.5;

        "ILambda.*"     0.5;

        Qr              0.1;//0.5;

        G               0.1;//0.5;

    }

}



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

I have tried to change the realaxationFactors and the solvers tolerances but I got nothing. It always crashes at the third iteration. I know it's something related to the use of Qr at the boundaries because I tried it without using Qr and it works fine (I haven't ran it till the end, but at least it doesn't crash at the third iteration).

I would like to get some help to be able to solve this problem. In case you need more information feel free to ask.

Thanks in advance,

Alex

Shouldn't it be like following?

Code:

        QrNbr           Qr;

        Qr              none;

Thanks for the hint Paritosh, but no, this is not what solves the problem. Remember I was talking about the external boundaries so no neighbour region can bring radiation to the domain.

Finally I could find out the proper set up for the externalWallHeatFluxTemperature boundary condition so that it can include Qr in the thermal balance. The correct use of this BC must be as follows:

Code:

QrName Qr;

Is it running for you now? I had a (possibly similar) issue with temperatureRadCoupledMixed BC, which requires a specific cell size next to the patch for convergence.

What is exactly your problem chriss85? Maybe it can be useful for me and other users to know about it to be able to face similar problems in the future.

@zfaraday: Did you make any changes to thermoType setup in thermophysicalProperties file for using Qr?

@vasava: I didn't make any change in the thermophysicalProperties file. Was I supposed to modify this file in order to use Qr? If so, what should I do within the file?

Quote:

Originally Posted by zfaraday (Post 517312)

@vasava: I didn't make any change in the thermophysicalProperties file. Was I supposed to modify this file in order to use Qr? If so, what should I do within the file?

I am not sure that is why I was asking. Since you could run your case without making modifications I guess you need not make any modifications.

Was there any solution to this issue? I have the same problems now ....

I have also tried the

Code:

QrName Qr;

as suggested by zfaraday, but this leads to

Code:

Qr none;

in later iterations. The correct syntax seems to be

Code:

Qr Qr;

Hello Tina,

Quote:

Originally Posted by TinaB (Post 547054)

Was there any solution to this issue? I have the same problems now ....

I have also tried the

Code:

QrName Qr;

as suggested by zfaraday, but this leads to

Code:

Qr none;

in later iterations. The correct syntax seems to be

Code:

Qr Qr;

Yes, you are right, the correct sintax is the one you suggested. I apologize for not reporting it when I found it out.

The problem is that the solver usually crashes when Qr is added to the patch's balance, I don't know why, I'm actually testing it's behavior to check if it should be reported as a bug or if I do something wrong. If you test it, please, tell us here if it works properly for you or, otherwise, the solver also crashes in your case.

Best regards,

Alex

Hello Alex,

thanks for the quick reply.

I am testing it, and also in my case the solver is crashing. I am currently investigating using the buoyantSimpleFoam solver (for my first tests, I used chtMultiRegionSimpleFoam, and the behaviour was the same).

I did the following:

1) copy the hotRadiationRoomFvDOM tutorial
2) modify the boundary condition in 0/T for the fixedWalls:

Code:

fixedWalls

    {

         type            externalWallHeatFluxTemperature;

        kappa           fluidThermo;

        kappaName       none;

        Qr             Qr;

        q               uniform 0;

        value           uniform 300;

    }

3) use the Allrun script.

First, the density started to get negative, and at T = 22 it crashed:

Code:

Time = 21



DILUPBiCG:  Solving for Ux, Initial residual = 0.0121552, Final residual = 3.62233e-05, No Iterations 1

DILUPBiCG:  Solving for Uy, Initial residual = 0.00880712, Final residual = 2.64742e-05, No Iterations 1

DILUPBiCG:  Solving for Uz, Initial residual = 0.00823497, Final residual = 4.85035e-05, No Iterations 1

DILUPBiCG:  Solving for h, Initial residual = 0.999913, Final residual = 0.00503933, No Iterations 1

DICPCG:  Solving for p_rgh, Initial residual = 0.999985, Final residual = 0.00628023, No Iterations 22

time step continuity errors : sum local = 8.14024, global = 3.09588e-15, cumulative = 3.57527e-15

rho max/min : 976.285 -86.0505

DILUPBiCG:  Solving for epsilon, Initial residual = 0.990103, Final residual = 0.00279567, No Iterations 1

DILUPBiCG:  Solving for k, Initial residual = 0.981696, Final residual = 0.022613, No Iterations 1

ExecutionTime = 100.32 s  ClockTime = 103 s



Time = 22



DILUPBiCG:  Solving for Ux, Initial residual = 0.155555, Final residual = 0.00101554, No Iterations 1

DILUPBiCG:  Solving for Uy, Initial residual = 0.175699, Final residual = 0.00091969, No Iterations 1

DILUPBiCG:  Solving for Uz, Initial residual = 0.110968, Final residual = 0.000687635, No Iterations 1

DILUPBiCG:  Solving for h, Initial residual = 0.0243916, Final residual = 2.11674e-15, No Iterations 1





--> 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

    in file /home/martina/OpenFOAM-Code/OpenFOAM-2.3.x/OpenFOAM-2.3.x/src/thermophysicalModels/specie/lnInclude/thermoI.H at line 76.



FOAM aborting

Can you reproduce this? Did I do something badly wrong?

I'd like to simulate cases with heat transfer by convection and radiation, and I would like to use externalWallHeatFluxTemperature in basically all of its available modes (fixed heat flux, fixed heat transfer coefficient), including radiation.

If, for a start, I want an adiabatic wall, the zeroGradient boundary condition as in the hotRadiationRoom tutorial does give zero convective heat flux at the wall, but it does not touch the radiative heat flux.

Can this be implemented in some other way, if externalWallHeatFluxTemperature fails to work with radiation?

posible disgnostic

Hi TIna,

Finally, today I've had time to test the BC for my simple case, described and uploaded here, and I got some answers about what is going wrong...

As you should know, externalWallHeatFlux BC works in two possible modes: constant heat flux(q) and constant heat transfer coefficient(h and Ta). In my case I only work with the second mode so I will describe what happens in the case I'm working on.

First of all, we must know what is the BC doing during the calulcation, so we must go to the source code. Here we can see that the values used for this BC are the following ones:

$refValue = \frac{hp*Ta}{hp- \frac{Qr}{Tp}}; refGrad = 0.0$
$f = \frac{hp-\frac{Qr}{Tp}}{hp-\frac{Qr}{Tp}+\frac{k}{\delta}}$

being

$hp= \frac{1}{\frac{1}{h}+\frac{\delta_e}{k_e}}$

the fluid's thermal conductivity and $\delta$ the cell center to patch' face center distance

with this values OpenFOAM calculates the temperature at the patch as

$Tp = f*refValue + (1-f)*[Ti + refGrad*\delta]$

In my case, I use very low values of cp for the air since I only want to compute radiative heat transfer, neglecting convection. Thus, the thermal conductivity for the air is close to 0 what gives a value for "f" of 1. Then, in my case the value of Tp is equal to "refValue". And here comes the problem!

The value of refValue becomes negative when $\frac{Qr}{Tp}>hp$ and

and, in my case, it happens in the 16th time step, so the solver crashes because of unphysical temperatures are reached at the patch.

I tried to modify the BC specification by changing the values mentioned above into the following ones

$refValue = Ta +\frac{Qr}{hp}; refGrad = 0.0$
$f = \frac{hp}{hp- \frac{k}{\delta}}$

that are supposed to compute the same balance at the patch so, in case that nothing crashed, both should give the same values. In this case, though, it also crashed at the 22th time step, so the matter is the same, Qr gives trouble when it is added to the balance because it gives negative values for the refValue term. This leads to negative (unphysical) temperatures at the patch so the solver stops showing an error.

@TinaB: In your case (constant heat flux) the values at the patch are computed with the following specification:

$refGrad= \frac{q + Qr}{k}; refValue = 0.0; f= 0.0$

It can be seen that when Qr gets negative values [and $\frac{|Qr|}{k}*\delta >Ti$ ]* the solver will also crash.

At this stage, I haven't found a solution that can be used in order to incorporate Qr to the balance at patches with this BC. Maybe something could be done by selecting accurately the initial values or the solver's tolerances and/or relaxation factors for some variables, but I don't know, just guessing... Maybe it is a bug that should be reportad to the bug tracker.

I will keep trying to solve it on my own when I have time, if someone els has faced the same problem than us and has been able to solve it somehow, please, let us know!

Best regards,

Alex

__________________________________________________ ____
Edit: added information between [ ]*

Hello Alex,

Thanks a lot for the detailed analysis. I was thinking along the same lines ...

This seems related to a bug report + bugfix concerning the turbulentTemperatureRadCoupledBaffleMixed: http://openfoam.org/mantisbt/view.php?id=1369

I am trying to work out some details. I was thinking of calculating the radiation temperature based on the incoming radiative heat flux Qin, and using this somehow in a mixed (or fixedValue) boundary condition and avoid using Qr as a gradient (or as part of the refValue) - though I don't know how (if) this will work out in detail.

I have the feeling that Qr is numerically unstable (maybe because it is calculated as the sum two large numbers Qin and Qem, but I did not find this in the code yesterday) - I have tried the hotRadiationRoom Tutorial, and set all the temperatures to 300 K such that basically no net radiation should occur. The calculation still diverged.

Best Regads,
Tina

possibly good news!

Hi Tina,

After a couple of trials I think this morning I have found a possible solution for our little problem! :D

The fact is that Qr seems to get unstable values at the first iterations that lead the solver to a crash. Well, since I also work with thermalBaffle1D BC in the simple test I shared in the link attached in my last post, and I already had to struggle a little with its source code, I decided to try using a relaxation factor for Qr like it is done in thermalBaffle1D.

I implemented it in the modified externalWallHeatFLuxTemperature BC I mentioned in my previous post and, after some tests I found out that the solver converged with a value for the relaxation factor of 0.1 (a very low value, but at least it made it converge). I don't know if the values I get after the solver converges are correct or not, but at least they look physically possible.

I still haven't tried if the solver converges by applying this modification in the original BC instead the one I modified. I hope it will too.

Now you can try to do what I suggested in order to use this BC to add Qr to the boundaries' balance and check if the results are nice or not. If you don't know how to do so contact me and I can send you the modified source code.

Best regards,

Alex

Hi Alex,

Congratulations, the idea is wonderful, and it works! Thanks for the Hint!

I tried it for the hotRadiationRoomFvDOM tutorial with the modifications I made, and now the energy balance includes radiation, and the total heat flux is (almost) the one that I specify in the boundary condition (q=0 in my case). The relaxation factor I used was 0.02; at 0.1 the calculation diverged, and at 0.05 the wall temperatures oscillated.

To check the results, I also calculated the case in Fluent, and the temperatures look very similar. I will perform further tests in other cases, and report my findings.

In the morning, I made some experiments myself: I tried to limit the gradient, such that if $\frac{Q_r}{\kappa}\delta +T_i \le T_{min}$ the BC switched to a fixedValue with $RefValue = T_{min}$ , and if $\frac{Q_r}{\kappa}\delta +T_i \ge T_{max}$ , I put $RefValue = T_{max}$ .

This did not diverge, but oscillated a lot and did not converge either.

I think that still a bug report should be filed, with the suggested fix.

Best Regards,
Tina

Hi Tina,

I'm glad you find it also useful! I still haven't tried it further enough to ensure its behavior now is correct, however, as per what you say, it seems it is a decent solution although I don't like the fact that the required values for the relaxation factor are that low...

I already reported the bug, and proposed the possible fix but I still haven't get an answer still.

Best regards,

Alex

Hello,

I am currently trying to model a solar volumetric receiver, i have conducted optical ray tracing analysis and have value for the heat flux at the aperture of the receiver. i have the data for the heat flux in 3 formats, grid form, matrix form, and table form. the heat flux represents a spot focus area (3cm radius) with energy density values in the formats mentioned. if you can shed some light on this issue it would be really helpful.

thanks