CFD Online Discussion Forums - chtMultiRegionFoam wallHeatFlux mismatch

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Post-Processing (https://www.cfd-online.com/Forums/openfoam-post-processing/)

- - chtMultiRegionFoam wallHeatFlux mismatch (https://www.cfd-online.com/Forums/openfoam-post-processing/221164-chtmultiregionfoam-wallheatflux-mismatch.html)

chtMultiRegionFoam wallHeatFlux mismatch

Dear CFD-ehtousiasts,

Currently I am working on a model that should capture the convective heat transfer between a coil, with hot water running through it, to ambient tankwater. With this model I try to capture the local and total heat transfer.

The first setup of the model includes a simple single tube surrounded by tankWater.

I am extracting heat flux values on both sides of my coil and in both directions (just for checking). These value nicely match.

Code:

wallHeatFlux1

{

    type        wallHeatFlux;

    libs        ("libfieldFunctionObjects.so");

    patches     ("coilWall_to_coilWater");

    region      coilWall;

    executeControl  writeTime;

    writeControl    writeTime;

}



wallHeatFlux2

{

    type        wallHeatFlux;

    libs        ("libfieldFunctionObjects.so");

    patches     ("coilWater_to_coilWall");

    region      coilWater;

    executeControl  writeTime;

    writeControl    writeTime;

}



wallHeatFlux3

{

    type        wallHeatFlux;

    libs        ("libfieldFunctionObjects.so");

    patches     ("coilWall_to_tankWater");

    region      coilWall;

    executeControl  writeTime;

    writeControl    writeTime;

}



wallHeatFlux4

{

    type        wallHeatFlux;

    libs        ("libfieldFunctionObjects.so");

    patches     ("tankWater_to_coilWall");

    region      tankWater;

    executeControl  writeTime;

    writeControl    writeTime;

}

Now, the problem:

If I compare the wall heat flux values from the wallHeatFlux utility, to the value computed from the flow rate and temperature at the inlet and outlet as:
Qdot=(mdot*rho*Cp*(inletT-outletT))

I find a significant mismatch between both values (see figure). BTW: I used patchAverage to compute my mean outletT.

Unfortunately, I could not find a solution in other threads on the wallHeatFlux utility.

Any ideas or tips would be much appreciated!

I ran the simulation with turbulence turned off. Surprisingly (or not) the wall heat flux and flowrate based heat flux to almost match perfectly (see figure).

Now, I am using an KOmegaSST turbulence model. In order to capture the heat transfer as good as possible I have meshed my system in such a way that y+ is smaller or equal to one everywhere and therefore not using any wall functions.

Any experience on this aspect what could cause the energy imbalance?

Your formula is only an approximation instead of temperature you'd need to look at the total Enthalpy between inlet and outlet. And you'd need to integrate instead of average everything.You would also need to include kinetic energy, although this part should probably nearly cancel out or be insignificant. Unless maybe the diameter changed drastically and the flow is quite fast.
Are any of your material or transport properties temperature dependent? Do the wall heat fluxes from the different regions match?
So is the coilWall_to_coilWater heatflux nearly identical to the coilWater_to_ coilWall flux?

Thanks for your help! As mentioned in the original post I checked the heat flux on each side and they do match. Regarding the averaging of the temperature, patchAverage uses areaAverage and therefore does a weighted average calculation. As far as I know this is equivalent to integration, right?

I understand your point about the total enthalpy. However I do not believe that the internal energy/velocity would change so significant that it results in the given energy mismatch. It is still worth checking. Best way would to this would be to check the enthalpy flux over the inlet and outlet? I will have to look into this.

I do use temperature dependent material properties, however I do not believe the problem can be found here, as the no-turbulence model gives good results with the same material properties.

BTW: tube is horizontal and has a uniform cross-section.

This is the full energy equation (not a single simplification made):

$0 = \nabla\cdot\left(\left(e+\frac{1}{2}\mathbf{v}\cdot\mathbf{v}\right)\rho\mathbf{v}\right)-\nabla\cdot\left(\lambda\nabla T\right)+\nabla\cdot\left(\mathbf{v}p\right)-\underbrace{\rho\mathbf{g}\cdot\mathbf{v}-\rho\dot{q}-\nabla\cdot\left(\boldsymbol{\tau}\cdot\mathbf{v}\right)}_{0}$
The last part of it assumed as zero since

gravitational acceleration is not active
no internal heat source is active inside the fluid
you are not adding viscous dissipation via fvOption. This is often negligible, but needs to be accounted for if you are simulating highly viscous substances.

With that said we are left with only these terms:
$0 = \nabla\cdot\left(\left(e+\frac{1}{2}\mathbf{v}\cdot\mathbf{v}\right)\rho\mathbf{v}\right)-\nabla\cdot\left(\lambda\nabla T\right)+\nabla\cdot\left(\mathbf{v}p\right)$

Now we apply the Gauss theorem:

$\int \nabla\cdot\mathbf{f}dV = \int \mathbf{n}\cdot\mathbf{f}dA$

Which yields:

$0 = \int \mathbf{n}\cdot\left(\left(e+\frac{1}{2}\mathbf{v}\cdot\mathbf{v}\right)\rho\mathbf{v}-\lambda\nabla T+\mathbf{v}p\right)\, dA$

This equals:

$0 = \int \mathbf{n}\cdot\left(H\rho\mathbf{v}-\lambda\nabla T\right)\, dA$

where

$H = h+\frac{1}{2}\mathbf{v}\cdot\mathbf{v}$

and the pv term was used to transform e to h. We now need to evaluate that term on the entire boundary of our domain. I'll split it into inlet/outlet and walls. With the difference being the assumption that nothing flows through the wall (n.v=0)

$\underbrace{\int_{wall} \mathbf{n}\cdot\left(\lambda\nabla T\right)\, dA_{wall}}_{\dot{Q}} = \int_{inout} \mathbf{n}\cdot\left(H\rho\mathbf{v}-\lambda\nabla T\right)\, dA_{inout}$

Under the assumption that no temperature gradient is present at inlet and outlet and no change in kinetic energy this equals your formula

$\dot{Q}= \int_{inout} \mathbf{n}\cdot\left(h\rho\mathbf{v}\right)\, dA_{inout}= \int_{inout} \mathbf{n}\cdot\left(c_p T \rho\mathbf{v}\right)\, dA_{inout}$

Now you can check if either of these assumptions is not the case for your simulation. And no you need to integrate since T is multiplied and hence weighted by v and cp if neither of these 3 is constant this is inaccurate.

Thank you very much! I now integrated the inlet and outlet field using patchMassFlowAverage

Code:

 massFlowAverageT

    {

    type patchMassFlowAverage;

    functionObjectLibs ( "libsimpleFunctionObjects.so" );

    fields ( T );

    patches

    (

        inlet

        outlet

    );

    factor 1.0;

    region coilWater;

    outputControlMode outputTime;

    outputInterval 1;

    writeStartTime no;

    verbose true;

    }

In this way you do account for the varying velocity field over the outflow patch together with temperature and indeed the wall heat flux and flow rate based Qdot match.

Naturally, if one wants to be really precise the variation of Cp over the patch should also be taken into account as you showed in the derivation. However, taking into account the minimal variation of Cp (and the hassle to get obtain that variable in post-processing) I left that out for now.

Again thanks a lot for your detailed elaboration!