heat not balanced in the chtMultiRegionSimpleFoam solver
 

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:
https://i.ibb.co/yWdCHDx/image.png
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:
https://i.ibb.co/bQVRqvQ/image.png
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!:o

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?

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%)...

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.

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

You have to calculate the temperature massflow averaged then it should be correct.

Hi, Adam!

Can you tell me how to calculate the temperature mass flow averaged?

I use something like this in my controlDict

Hi, Adam!

I tried the weightedAverage and now...I'm confused...: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.34e-4*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...

Are youre transport properties constant?

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:

Have you tried the function object for the energy summation (weightedSum)? This should be the energy flow in J/s through the inlet/outlet.

Yes I tried it, the phi weighted Sum of h, and the result was on #5.
The energy change is not balanced.

Is your result mesh independent? If not try a finer mesh and see if you still have this behavior.

And what are the residuals of p (p_rgh?) and U? The initial residual of the last iteration should be below 1e-4 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?

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.470719e-06, No Iterations 1
DILUPBiCGStab: Solving for Uy, Initial residual = 0.0007633498, Final residual = 1.238344e-06, No Iterations 1
DILUPBiCGStab: Solving for Uz, Initial residual = 0.001816722, Final residual = 3.271993e-06, No Iterations 1

DILUPBiCGStab: Solving for h, Initial residual = 0.0009938677, Final residual = 5.421818e-06, No Iterations 1

DICPCG: Solving for p_rgh, Initial residual = 0.004760268, Final residual = 4.204527e-05, No Iterations 14
DICPCG: Solving for p_rgh, Initial residual = 8.245556e-05, Final residual = 7.71274e-07, No Iterations 501
DICPCG: Solving for p_rgh, Initial residual = 4.175771e-06, Final residual = 8.664596e-08, No Iterations 17
DICPCG: Solving for p_rgh, Initial residual = 3.342495e-07, Final residual = 7.87855e-08, No Iterations 2

time step continuity errors : sum local = 0.0002149379, global = 1.883727e-05, 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...

The only idea I have: what are the temperature boundary conditions of the walls of the collant region? I guess all should be zeroGradient.

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

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!

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 1e-3 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 1e-6 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~