CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Running, Solving & CFD

Problems adding Qr field in externalWallHeatFluxTemperature BC

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree4Likes
  • 2 Post By zfaraday
  • 1 Post By zfaraday
  • 1 Post By zfaraday

Reply
 
LinkBack Thread Tools Display Modes
Old   October 29, 2014, 16:04
Exclamation Problems adding Qr field in externalWallHeatFluxTemperature BC
  #1
Senior Member
 
Alex
Join Date: Oct 2013
Posts: 267
Rep Power: 12
zfaraday will become famous soon enough
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
Attached Images
File Type: jpg geometry.jpg (35.3 KB, 22 views)
Attached Files
File Type: txt log.chtMultiRegionSimpleFoam.txt (71.8 KB, 6 views)
scintilla and flowAlways like this.
__________________
I'm newbie in OpenFOAM's world and not an English-speaking, so if I make any mistake a correction will be welcome!
zfaraday is online now   Reply With Quote

Old   November 4, 2014, 06:57
Default
  #2
Senior Member
 
Paritosh Vasava
Join Date: Oct 2012
Location: Lappeenranta, Finland
Posts: 544
Rep Power: 12
vasava will become famous soon enough
Shouldn't it be like following?

Code:
        QrNbr           Qr;
        Qr              none;
vasava is offline   Reply With Quote

Old   November 4, 2014, 07:19
Default
  #3
Senior Member
 
Alex
Join Date: Oct 2013
Posts: 267
Rep Power: 12
zfaraday will become famous soon enough
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;
scintilla likes this.
__________________
I'm newbie in OpenFOAM's world and not an English-speaking, so if I make any mistake a correction will be welcome!
zfaraday is online now   Reply With Quote

Old   November 4, 2014, 07:23
Default
  #4
Senior Member
 
Join Date: Oct 2013
Posts: 237
Rep Power: 4
chriss85 is on a distinguished road
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.
chriss85 is offline   Reply With Quote

Old   November 4, 2014, 07:38
Default
  #5
Senior Member
 
Alex
Join Date: Oct 2013
Posts: 267
Rep Power: 12
zfaraday will become famous soon enough
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.
__________________
I'm newbie in OpenFOAM's world and not an English-speaking, so if I make any mistake a correction will be welcome!
zfaraday is online now   Reply With Quote

Old   November 4, 2014, 08:51
Default
  #6
Senior Member
 
Paritosh Vasava
Join Date: Oct 2012
Location: Lappeenranta, Finland
Posts: 544
Rep Power: 12
vasava will become famous soon enough
@zfaraday: Did you make any changes to thermoType setup in thermophysicalProperties file for using Qr?
vasava is offline   Reply With Quote

Old   November 4, 2014, 09:02
Default
  #7
Senior Member
 
Alex
Join Date: Oct 2013
Posts: 267
Rep Power: 12
zfaraday will become famous soon enough
@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'm newbie in OpenFOAM's world and not an English-speaking, so if I make any mistake a correction will be welcome!
zfaraday is online now   Reply With Quote

Old   November 4, 2014, 09:21
Default
  #8
Senior Member
 
Paritosh Vasava
Join Date: Oct 2012
Location: Lappeenranta, Finland
Posts: 544
Rep Power: 12
vasava will become famous soon enough
Quote:
Originally Posted by zfaraday View Post
@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.
vasava is offline   Reply With Quote

Old   November 4, 2014, 09:59
Default
  #9
Senior Member
 
Join Date: Oct 2013
Posts: 237
Rep Power: 4
chriss85 is on a distinguished road
See Heat transfer convergence problem and Understanding temperature coupling BCs.
chriss85 is offline   Reply With Quote

Old   May 20, 2015, 10:44
Default
  #10
New Member
 
Martina Blank
Join Date: Feb 2014
Posts: 7
Rep Power: 3
TinaB is on a distinguished road
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;
.
TinaB is offline   Reply With Quote

Old   May 20, 2015, 10:53
Default
  #11
Senior Member
 
Alex
Join Date: Oct 2013
Posts: 267
Rep Power: 12
zfaraday will become famous soon enough
Hello Tina,
Quote:
Originally Posted by TinaB View Post
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
__________________
I'm newbie in OpenFOAM's world and not an English-speaking, so if I make any mistake a correction will be welcome!
zfaraday is online now   Reply With Quote

Old   May 20, 2015, 11:27
Default
  #12
New Member
 
Martina Blank
Join Date: Feb 2014
Posts: 7
Rep Power: 3
TinaB is on a distinguished road
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?
TinaB is offline   Reply With Quote

Old   May 21, 2015, 15:28
Exclamation posible disgnostic
  #13
Senior Member
 
Alex
Join Date: Oct 2013
Posts: 267
Rep Power: 12
zfaraday will become famous soon enough
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}}
k 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 Qr>0 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 [ ]*
flowAlways likes this.
__________________
I'm newbie in OpenFOAM's world and not an English-speaking, so if I make any mistake a correction will be welcome!

Last edited by zfaraday; May 21, 2015 at 17:19. Reason: added information missing
zfaraday is online now   Reply With Quote

Old   May 22, 2015, 02:58
Default
  #14
New Member
 
Martina Blank
Join Date: Feb 2014
Posts: 7
Rep Power: 3
TinaB is on a distinguished road
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
TinaB is offline   Reply With Quote

Old   May 22, 2015, 07:59
Cool possibly good news!
  #15
Senior Member
 
Alex
Join Date: Oct 2013
Posts: 267
Rep Power: 12
zfaraday will become famous soon enough
Hi Tina,

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

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
__________________
I'm newbie in OpenFOAM's world and not an English-speaking, so if I make any mistake a correction will be welcome!
zfaraday is online now   Reply With Quote

Old   May 22, 2015, 10:36
Default
  #16
New Member
 
Martina Blank
Join Date: Feb 2014
Posts: 7
Rep Power: 3
TinaB is on a distinguished road
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
TinaB is offline   Reply With Quote

Old   May 22, 2015, 12:58
Default
  #17
Senior Member
 
Alex
Join Date: Oct 2013
Posts: 267
Rep Power: 12
zfaraday will become famous soon enough
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
__________________
I'm newbie in OpenFOAM's world and not an English-speaking, so if I make any mistake a correction will be welcome!
zfaraday is online now   Reply With Quote

Reply

Tags
baffle, chtmultiregion, frozenflow, radiation, vacuum

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Foam::error::PrintStack almir OpenFOAM Running, Solving & CFD 51 June 28, 2015 16:36
dynamicTopoFVMesh and pointDisplacement RandomUser OpenFOAM Meshing & Mesh Conversion 5 March 2, 2015 18:19
Is good initial guess field is neccessary ? mmkr825 OpenFOAM 5 October 17, 2012 12:58
adding temperature field to icofoam houkensjtu OpenFOAM 5 September 26, 2012 22:20
Adding a user application class Rasmus Gjesing (Gjesing) OpenFOAM Pre-Processing 57 February 3, 2010 04:45


All times are GMT -4. The time now is 13:51.