
[Sponsors] 
heat not balanced in the chtMultiRegionSimpleFoam solver 

LinkBack  Thread Tools  Search this Thread  Display Modes 
June 30, 2019, 21:59 
heat not balanced in the chtMultiRegionSimpleFoam solver

#1 
Member
Join Date: May 2013
Posts: 31
Rep Power: 9 
Hi, foamers,
I met a difficult problem that in my calculation of a single shell and tube reactor the heat is not balanced Briefly, I used chtMultiRegionSimpleFoam. The mesh had 3 regions  gas region, reactor tube, and coolant region. For all the regions the schemes were set as limitedLinear 1 for div and linear corrected for laplacian. The gas at the inlet and outlet had the same temperature, thus I thought that all the reaction heat was passed and brought out by the coolant. The reaction heat release was calculated in the solver by fvc::domainIntegrate(Qdot). The heat brought by the coolant was calculated using paraview by Q=C_p*m*delta_T using the temperature difference between the coolant inlet and outlet. The gas region and reactor tube, reactor tube and coolant region were coupled by compressible::turbulentTemperatureRadCoupledMixed. However, the problem happens  the released heat by the reaction was not balanced with the heat brought out by the coolant! Here I post this maybe you can give me some advises or direct my problem... Next is the detail. First I used a mesh like this: the red part is coolant region, green is the reactor tube that only heat transfer is calculated, white is gas region. the calculation error is like this: divSchemes of coolant, reaction heat, heat passed by the turbulentTemperatureRadCoupledMixed, heat brought out by coolant, error bounded Gauss limitedLinear 1, 975J/s, 975J/s , 939J/s, 3.69% bounded Gauss linearUpwind grad(U), 975J/s, 975J/s , 947J/s, 2.87% Gauss upwind, 975J/s, 975J/s , 963J/s, 1.13% Because the reaction heat the heat passed by the turbulentTemperatureRadCoupledMixed between the reactor tube and coolant, I thought that the error comes from the coolant calculation... And another is that the upwind of div gives a more accurate result... I hope you can figure out why this happens. So, anyway I think the error may related to the coolant mesh, I used a mesh like then: It does not have the coolant in and out pipe, the coolant directly flows through the left and right surface. the calculation error is like this: divSchemes of coolant, reaction heat, heat passed by the turbulentTemperatureRadCoupledMixed, heat brought out by coolant, error bounded Gauss linear, 962J/s, 962J/s , 958J/s, 0.42% Gauss upwind, 962J/s, 962J/s , 960J/s, 0.21% Wow, the error is quite small. However I cannot use a mesh like this because the flow pattern is different with and without the coolant in,out pipe. Above calculation told me that the mesh with coolant in,out pipe maybe not good enough, so I did this: using the mesh with coolant in,out pipe but remove the gas region(white part), and set the inner surface of the reactor tube as externalWallHeatFluxTemperature, with power of Q 964 J/s. the calculation error is like this: divSchemes of coolant, externalWallHeatFluxTemperature, heat passed by the turbulentTemperatureRadCoupledMixed, heat brought out by coolant, error bounded Gauss limitedLinear 1, 964J/s, 964J/s , 968J/s, 0.44% ......It gave me a small error! but why??? I also made the mesh with snappyHexMesh without the final step of "snap"(which means a staggered hex structured mesh), but the result is the same, if I include the gas region, then the coolant heat will not balance, and the error becomes larger when I use limitedLinear scheme. So, as conclusion of my problem, if I calculate the mesh include the gas region, 1. reaction heat > 2. surface between gas and tube(turbulentTemperatureRadCoupledMixed) > 3. surface between tube and coolant(turbulentTemperatureRadCoupledMixed) > 4. coolant brought out until 3, the heat was balanced, the heat released by the reaction and passed by the coupled boundary was the same, only 4 was not balanced. If I remove the gas region, 2. surface between gas and tube(externalWallHeatFluxTemperature, given power Q) > 3. surface between tube and coolant(turbulentTemperatureRadCoupledMixed) > 4. coolant brought out then all the process had the heat balanced. Sorry for the long text, but I really hope that you can figure out my problem or give me some advises. Thank you! 

July 3, 2019, 07:01 

#2 
Senior Member
Joachim Herb
Join Date: Sep 2010
Posts: 567
Rep Power: 16 
Have you tried to use the wallHeatFlux function object to calculate the heat flux at the inner and outer surface of your reactor tube. What are the results? Are they the same on both surfaces and on the inner and outer side of each surface, meaning if you do the calculation for the different regions?


July 3, 2019, 22:20 

#3 
Member
Join Date: May 2013
Posts: 31
Rep Power: 9 
Hi, Herb!
Yes, the heat flux is the same between the inner and outer side of the surface. First I turned on the debug switch for the "compressible::turbulentTemperatureRadCoupledMixed " boundary so it outputs the heat transfer rate during the calcualtion. It said like: coolant_region:coolant_region_to_solid_pipe_region :T < solid_pipe_region:solid_pipe_region_to_coolant_reg ion:T : heat transfer rate:964.9667 walltemperature min:473.149 max:494.9711 avg:475.3582 for the coolant and said like solid_pipe_region:solid_pipe_region_to_coolant_reg ion:T < coolant_region:coolant_region_to_solid_pipe_region :T : heat transfer rate:964.9698 walltemperature min:473.149 max:494.97 avg:475.3582 for the tube region. The heat flux is the same. I also did your suggestion to use wallHeatFlux function object, and it said like this wallHeatFlux wallHeatFlux_coolantToTube write: writing field wallHeatFlux min/max/integ(coolant_wall) = 0, 0, 0 min/max/integ(coolant_region_to_solid_pipe_region) = 653.3906, 48224.13, 964.9723 wallHeatFlux wallHeatFlux_TubeToCoolant write: writing field wallHeatFlux min/max/integ(metal_wall) = 0, 0, 0 min/max/integ(solid_pipe_region_to_gas_region) = 1.193805, 57428.32, 965.064 min/max/integ(solid_pipe_region_to_coolant_region) = 48158.81, 653.3989, 964.9513 still the heat flux is the same... At this time, the heat brought out by the coolant outlet is calculated as 900 J/s, quite a large error(6.7%)... 

July 4, 2019, 07:00 

#4 
Senior Member
Joachim Herb
Join Date: Sep 2010
Posts: 567
Rep Power: 16 
You are calculating the heat flux through the outlet in paraview? Are you sure, the velocity values and the temperatures are correct? Are you using the values of the patches (have you loaded them explicitly) or is paraview using some interpolation of cell values?
Have you tried function objects to do the heat flux calculation in OpenFOAM directly? You could e.g. sum the field phi over the in/outlet to get the mass and energy flux calculated by OpenFOAM. Code:
massInlet { type surfaceFieldValue; libs ("libfieldFunctionObjects.so"); writeControl timeStep; writeInterval 1; log yes; writeTotalArea no; writeFields no; regionType patch; name inlet; operation sum; fields ( phi ); } energyInlet { type surfaceFieldValue; libs ("libfieldFunctionObjects.so"); writeControl timeStep; writeInterval 1; log yes; writeTotalArea no; writeFields no; regionType patch; name inlet; operation weightedSum; weightField phi; fields ( h ); } 

July 5, 2019, 05:27 

#5 
Member
Join Date: May 2013
Posts: 31
Rep Power: 9 
Hi, Herb!
This is how I calculate the energy brought out by coolant using paraview: 1. load the inlet boundary separately 2. apply the integral filter to the boundary, then I can get the area, averaged temperature, density, and velocity 3. repeat 1 and 2 for the outlet boundary 4. I can calculated the mass flux and temperature difference use these values, and finally the energy change by Q=C*m*deltaT I tried your suggestion of "surfaceFieldValue" function, Its results is like: coolant inlet mass flux = 0.64558 kg/s coolant outlet mass flux = 0.64558 kg/s coolant inlet energy flux = 664672.43 J/s (does this unit right?) coolant outlet energy flux = 665599.7428 J/s the energy difference is about 927.31 J/s (the error is about 3.8%, compared to the right value of 965 J/s) It is still not balanced to the heat transferred into the coolant through the coupled patch... 

July 5, 2019, 05:40 

#6 
New Member
Adam
Join Date: Jan 2019
Posts: 21
Rep Power: 3 
You have to calculate the temperature massflow averaged then it should be correct.


July 5, 2019, 05:53 

#7 
Member
Join Date: May 2013
Posts: 31
Rep Power: 9 
Hi, Adam!
Can you tell me how to calculate the temperature mass flow averaged? 

July 5, 2019, 06:14 

#8 
New Member
Adam
Join Date: Jan 2019
Posts: 21
Rep Power: 3 
Code:
INLET_TEMP { type surfaceFieldValue; libs ("libfieldFunctionObjects.so"); writeControl writeTime; log yes; writeFields no; regionType patch; name INLET; operation weightedAverage; weightField phi; surfaceFormat foam; fields ( T ); } I use something like this in my controlDict 

July 7, 2019, 22:04 

#9 
Member
Join Date: May 2013
Posts: 31
Rep Power: 9 
Hi, Adam!
I tried the weightedAverage and now...I'm confused... On both inlet and outlet boundary of coolant, I used sum of phi to calculate the mass flux, and it gave mass flux of inlet = outlet = 0.64558 kg/s use phi weightedAverage to calculate the temperature of each boundary and it gave T of inlet = 473.15 K T of outlet = 473.8074 K thus the delta_T is 0.6574 K use phi weightedAverage to calculate the velocity of each boundary and it gave U of inlet = 2.1 m/s U of outlet = 2.183 m/s the inlet and outlet area are the same. So, has above data I can calculate the heat flux by Q=C_p*m*delta_T=2176*0.64558*0.6574=923.5 J/s the error is still about 4.2%( right value 964 J/s ) but, If I calculate the heat flux by Q=C_p*Area*U*rho*delta_T=2176*3.34e4*2.183*920*0.6574=960 J/s then it almost the right value! From instinct I think the first calculation is right because it calculates the mass flux by openfoam. the second one's velocity is not the same at inlet and outlet... 

July 8, 2019, 03:09 

#10 
New Member
Adam
Join Date: Jan 2019
Posts: 21
Rep Power: 3 
Are youre transport properties constant?


July 8, 2019, 06:05 

#11 
Member
Join Date: May 2013
Posts: 31
Rep Power: 9 
yes, they are the same.
At first, I did not realize that openfoam can calculate the flux for me so I set them all constant that I can use paraview to process them. here the thermophysicalProperties for coolant: Code:
thermoType { type heRhoThermo; mixture pureMixture; transport const; thermo hConst; equationOfState rhoConst; specie specie; energy sensibleEnthalpy; } 

July 8, 2019, 09:15 

#12 
Senior Member
Joachim Herb
Join Date: Sep 2010
Posts: 567
Rep Power: 16 
Have you tried the function object for the energy summation (weightedSum)? This should be the energy flow in J/s through the inlet/outlet.


July 8, 2019, 21:25 

#13 
Member
Join Date: May 2013
Posts: 31
Rep Power: 9 
Yes I tried it, the phi weighted Sum of h, and the result was on #5.
The energy change is not balanced. 

July 9, 2019, 03:18 

#14 
New Member
Adam
Join Date: Jan 2019
Posts: 21
Rep Power: 3 
Is your result mesh independent? If not try a finer mesh and see if you still have this behavior.


July 9, 2019, 18:14 

#15 
Senior Member
Joachim Herb
Join Date: Sep 2010
Posts: 567
Rep Power: 16 
And what are the residuals of p (p_rgh?) and U? The initial residual of the last iteration should be below 1e4 or even smaller. (Normally the residuals of h are smaller anyway.)
Is your simulation using a turbulence model? Have you tried laminar? Have you tried a transient simulation? 

August 9, 2019, 03:15 

#16 
Member
Join Date: May 2013
Posts: 31
Rep Power: 9 
Hi! sorry for the late reply.
I made a finer mesh of the coolant region and tested it. the original one has about 170000 cells and this new one has about 370000 cells. the result is the same, the heat is not balanced... both RAS and laminar were tried, the heat is not balanced... below is the residuals of the last step: DILUPBiCGStab: Solving for Ux, Initial residual = 0.002163327, Final residual = 3.470719e06, No Iterations 1 DILUPBiCGStab: Solving for Uy, Initial residual = 0.0007633498, Final residual = 1.238344e06, No Iterations 1 DILUPBiCGStab: Solving for Uz, Initial residual = 0.001816722, Final residual = 3.271993e06, No Iterations 1 DILUPBiCGStab: Solving for h, Initial residual = 0.0009938677, Final residual = 5.421818e06, No Iterations 1 DICPCG: Solving for p_rgh, Initial residual = 0.004760268, Final residual = 4.204527e05, No Iterations 14 DICPCG: Solving for p_rgh, Initial residual = 8.245556e05, Final residual = 7.71274e07, No Iterations 501 DICPCG: Solving for p_rgh, Initial residual = 4.175771e06, Final residual = 8.664596e08, No Iterations 17 DICPCG: Solving for p_rgh, Initial residual = 3.342495e07, Final residual = 7.87855e08, No Iterations 2 time step continuity errors : sum local = 0.0002149379, global = 1.883727e05, cumulative = 0.1843702 although I made a transient solver, I did not use it because it needs a long time to calculate to reach the steady state... 

August 15, 2019, 07:44 

#17 
Senior Member
Joachim Herb
Join Date: Sep 2010
Posts: 567
Rep Power: 16 
The only idea I have: what are the temperature boundary conditions of the walls of the collant region? I guess all should be zeroGradient.


August 17, 2019, 09:41 

#18 
Senior Member
Joachim Herb
Join Date: Sep 2010
Posts: 567
Rep Power: 16 
And another "debugging" option: What happens, if you turn of the solid and the gas fluid zone in chtMultiRegionFoam?
(Just another guess: do you use "dpdt off" for the coolant zone?) 

August 20, 2019, 02:49 

#19 
New Member
shach
Join Date: Apr 2019
Posts: 26
Rep Power: 3 
Hello,
Have you solved this problem? I am facing a similar problem. The heat flux is balanced when the material properties are constant. But it not balanced when the material properties, c_p, kappa, rho .etc are variables respect to temperature. Do you have any suggestion? Thanks! 

September 26, 2019, 04:25 

#20 
Member
Join Date: May 2013
Posts: 31
Rep Power: 9 
Hi, shach
Finally I "solved" the problem. I found that the case is not converged when the heat is not balanced. I think this maybe related to the mesh I used. The initial residual for h is always about 1e3 level even 10 thousand step is calculated. After I made a new mesh, and set the divSchemes for the coolant region to linearUpwind, the initial residual for h then became smaller than 1e6 soon, and the heat calculated by the phi weighted Sum of h at inlet/outlet is finally balanced. And I also found a very stable set of divSchemes: div(phi,U) bounded Gauss linearUpwindV grad(U); div(phi,K) bounded Gauss linearUpwind grad(U); div(phi,h) Gauss upwind; div(phi,epsilon) bounded Gauss linearUpwind grad(U); div(phi,k) bounded Gauss linearUpwind grad(U); I found this set is very stable. In all my calculations, I will first try the limitedLinear div scheme, if it won't converge, I will use this set and it always give me a converged results! Hope this helps you~ 

Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Which solver for combined cfd and solid state heat transfer?  heliosoph  Main CFD Forum  3  November 11, 2015 15:56 
Separate Convective and Radiative heat transfer in CFD post using fluent as a solver  Lemanes  FLUENT  1  July 6, 2015 10:31 
Heat and Moisture solver  Neumann time varing BC  Ricardo  OpenFOAM Programming & Development  1  August 6, 2013 14:13 
Error finding variable "THERMX"  sunilpatil  CFX  8  April 26, 2013 07:00 
solver radiative and convective heat fluxes in CFDPost  user  CFX  3  August 3, 2010 08:21 