
[Sponsors] 
October 29, 2014, 16:04 
Problems adding Qr field in externalWallHeatFluxTemperature BC

#1 
Senior Member
Alex
Join Date: Oct 2013
Posts: 278
Rep Power: 12 
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; } 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; } // ************************************************************************* // 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 1e5;//1e7; relTol 0;//0.01; } Ii { solver GAMG; tolerance 1e4; relTol 0; smoother symGaussSeidel; cacheAgglomeration true; nCellsInCoarsestLevel 10; agglomerator faceAreaPair; mergeLevels 1; maxIter 5; nPreSweeps 0; nPostSweeps 1; } /* { $h; tolerance 1e4;//1e7; relTol 0.01; } */ } SIMPLE { momentumPredictor off; nNonOrthogonalCorrectors 0; pRefCell 0; pRefValue 1e2;//100000; rhoMin rhoMin [1 3 0 0 0] 1e10; rhoMax rhoMax [1 3 0 0 0] 20; frozenFlow true; residualControl { /* h 1e3; G 1e3; "ILambda.*" 1e3; */ } } relaxationFactors { fields { } equations { h 0.3;//0.5; "ILambda.*" 0.5; Qr 0.1;//0.5; G 0.1;//0.5; } } // ************************************************************************* // 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
__________________
I'm newbie in OpenFOAM's world and not an Englishspeaking, so if I make any mistake a correction will be welcome! 

November 4, 2014, 06:57 

#2 
Senior Member
Paritosh Vasava
Join Date: Oct 2012
Location: Lappeenranta, Finland
Posts: 545
Rep Power: 12 
Shouldn't it be like following?
Code:
QrNbr Qr; Qr none; 

November 4, 2014, 07:19 

#3 
Senior Member
Alex
Join Date: Oct 2013
Posts: 278
Rep Power: 12 
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;
__________________
I'm newbie in OpenFOAM's world and not an Englishspeaking, so if I make any mistake a correction will be welcome! 

November 4, 2014, 07:23 

#4 
Senior Member
Join Date: Oct 2013
Posts: 270
Rep Power: 4 
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.


November 4, 2014, 07:38 

#5 
Senior Member
Alex
Join Date: Oct 2013
Posts: 278
Rep Power: 12 
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 Englishspeaking, so if I make any mistake a correction will be welcome! 

November 4, 2014, 08:51 

#6 
Senior Member
Paritosh Vasava
Join Date: Oct 2012
Location: Lappeenranta, Finland
Posts: 545
Rep Power: 12 
@zfaraday: Did you make any changes to thermoType setup in thermophysicalProperties file for using Qr?


November 4, 2014, 09:02 

#7 
Senior Member
Alex
Join Date: Oct 2013
Posts: 278
Rep Power: 12 
@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 Englishspeaking, so if I make any mistake a correction will be welcome! 

November 4, 2014, 09:21 

#8 
Senior Member
Paritosh Vasava
Join Date: Oct 2012
Location: Lappeenranta, Finland
Posts: 545
Rep Power: 12 
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.


November 4, 2014, 09:59 

#9 
Senior Member
Join Date: Oct 2013
Posts: 270
Rep Power: 4 

May 20, 2015, 10:44 

#10 
New Member
Martina Blank
Join Date: Feb 2014
Posts: 7
Rep Power: 3 
Was there any solution to this issue? I have the same problems now ....
I have also tried the Code:
QrName Qr; Code:
Qr none; Code:
Qr Qr; 

May 20, 2015, 10:53 

#11  
Senior Member
Alex
Join Date: Oct 2013
Posts: 278
Rep Power: 12 
Hello Tina,
Quote:
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 Englishspeaking, so if I make any mistake a correction will be welcome! 

May 20, 2015, 11:27 

#12 
New Member
Martina Blank
Join Date: Feb 2014
Posts: 7
Rep Power: 3 
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; } 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.62233e05, No Iterations 1 DILUPBiCG: Solving for Uy, Initial residual = 0.00880712, Final residual = 2.64742e05, No Iterations 1 DILUPBiCG: Solving for Uz, Initial residual = 0.00823497, Final residual = 4.85035e05, 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.09588e15, cumulative = 3.57527e15 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.11674e15, 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/OpenFOAMCode/OpenFOAM2.3.x/OpenFOAM2.3.x/src/thermophysicalModels/specie/lnInclude/thermoI.H at line 76. FOAM aborting 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? 

May 21, 2015, 15:28 
posible disgnostic

#13 
Senior Member
Alex
Join Date: Oct 2013
Posts: 278
Rep Power: 12 
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: being the fluid's thermal conductivity and the cell center to patch' face center distance with this values OpenFOAM calculates the temperature at the patch as 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 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 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: It can be seen that when Qr gets negative values [and ]* 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 [ ]*
__________________
I'm newbie in OpenFOAM's world and not an Englishspeaking, 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 

May 22, 2015, 02:58 

#14 
New Member
Martina Blank
Join Date: Feb 2014
Posts: 7
Rep Power: 3 
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 

May 22, 2015, 07:59 
possibly good news!

#15 
Senior Member
Alex
Join Date: Oct 2013
Posts: 278
Rep Power: 12 
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 Englishspeaking, so if I make any mistake a correction will be welcome! 

May 22, 2015, 10:36 

#16 
New Member
Martina Blank
Join Date: Feb 2014
Posts: 7
Rep Power: 3 
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 the BC switched to a fixedValue with , and if , I put . 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 

May 22, 2015, 12:58 

#17 
Senior Member
Alex
Join Date: Oct 2013
Posts: 278
Rep Power: 12 
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 Englishspeaking, so if I make any mistake a correction will be welcome! 

Tags 
baffle, chtmultiregion, frozenflow, radiation, vacuum 
Thread Tools  
Display Modes  


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 PreProcessing  57  February 3, 2010 04:45 