CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM

heat not balanced in the chtMultiRegionSimpleFoam solver

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

Like Tree1Likes
  • 1 Post By boundary93

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   June 30, 2019, 21:59
Question heat not balanced in the chtMultiRegionSimpleFoam solver
  #1
New Member
 
Join Date: May 2013
Posts: 13
Rep Power: 8
carye is on a distinguished road
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!
carye is offline   Reply With Quote

Old   July 3, 2019, 07:01
Default
  #2
Senior Member
 
Joachim Herb
Join Date: Sep 2010
Posts: 512
Rep Power: 14
jherb is on a distinguished road
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?
jherb is offline   Reply With Quote

Old   July 3, 2019, 22:20
Default
  #3
New Member
 
Join Date: May 2013
Posts: 13
Rep Power: 8
carye is on a distinguished road
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%)...
carye is offline   Reply With Quote

Old   July 4, 2019, 07:00
Default
  #4
Senior Member
 
Joachim Herb
Join Date: Sep 2010
Posts: 512
Rep Power: 14
jherb is on a distinguished road
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
        );
    }
jherb is offline   Reply With Quote

Old   July 5, 2019, 05:27
Default
  #5
New Member
 
Join Date: May 2013
Posts: 13
Rep Power: 8
carye is on a distinguished road
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...
carye is offline   Reply With Quote

Old   July 5, 2019, 05:40
Default
  #6
New Member
 
Adam
Join Date: Jan 2019
Posts: 9
Rep Power: 2
boundary93 is on a distinguished road
You have to calculate the temperature massflow averaged then it should be correct.
boundary93 is offline   Reply With Quote

Old   July 5, 2019, 05:53
Default
  #7
New Member
 
Join Date: May 2013
Posts: 13
Rep Power: 8
carye is on a distinguished road
Hi, Adam!

Can you tell me how to calculate the temperature mass flow averaged?
carye is offline   Reply With Quote

Old   July 5, 2019, 06:14
Default
  #8
New Member
 
Adam
Join Date: Jan 2019
Posts: 9
Rep Power: 2
boundary93 is on a distinguished road
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
Tobi likes this.
boundary93 is offline   Reply With Quote

Old   July 7, 2019, 22:04
Default
  #9
New Member
 
Join Date: May 2013
Posts: 13
Rep Power: 8
carye is on a distinguished road
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.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...
carye is offline   Reply With Quote

Old   July 8, 2019, 03:09
Default
  #10
New Member
 
Adam
Join Date: Jan 2019
Posts: 9
Rep Power: 2
boundary93 is on a distinguished road
Are youre transport properties constant?
boundary93 is offline   Reply With Quote

Old   July 8, 2019, 06:05
Default
  #11
New Member
 
Join Date: May 2013
Posts: 13
Rep Power: 8
carye is on a distinguished road
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;
}
carye is offline   Reply With Quote

Old   July 8, 2019, 09:15
Default
  #12
Senior Member
 
Joachim Herb
Join Date: Sep 2010
Posts: 512
Rep Power: 14
jherb is on a distinguished road
Have you tried the function object for the energy summation (weightedSum)? This should be the energy flow in J/s through the inlet/outlet.
jherb is offline   Reply With Quote

Old   July 8, 2019, 21:25
Default
  #13
New Member
 
Join Date: May 2013
Posts: 13
Rep Power: 8
carye is on a distinguished road
Yes I tried it, the phi weighted Sum of h, and the result was on #5.
The energy change is not balanced.
carye is offline   Reply With Quote

Old   July 9, 2019, 03:18
Default
  #14
New Member
 
Adam
Join Date: Jan 2019
Posts: 9
Rep Power: 2
boundary93 is on a distinguished road
Is your result mesh independent? If not try a finer mesh and see if you still have this behavior.
boundary93 is offline   Reply With Quote

Old   July 9, 2019, 18:14
Default
  #15
Senior Member
 
Joachim Herb
Join Date: Sep 2010
Posts: 512
Rep Power: 14
jherb is on a distinguished road
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?
jherb is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
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
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 CFD-Post user CFX 3 August 3, 2010 08:21


All times are GMT -4. The time now is 12:29.