Heat imbalance at multiregional boundary in chtMultiRegionFoam

Summary

When specific heat has a temperature dependence, heat flux imbalance occurs at the multiregional boundary.
I confirmed this phenomenon using wall heat flux in Field function objects.

※This imbalance occurs in not only compressible::turbulentTemperatureRadCoupledMixed but also in * externalWallHeatFluxTemperature*

Steps to reproduce

1. Download the attached file(one is constant Cp, another is temperature dependent Cp)
2. run chtMultiReigonFoam in each directory
3. compare postprocessing/l/wallHeatFlux1/0.0/wallHeatFlux.dat and postprocessing/r/wallHeatFlux1/0.0/wallHeatFlux.dat

Example case

imbalancedHeatFluxWhenTemperatureDependentCp.tar.gz

balancedHeatFluxWhenTemperatureIndependentCp.tar.gz

What is the current bug behaviour?

Heat flux imbalance occurs at the multiregional boundary.

What is the expected correct behavior?

Heat flux is balanced at the multiregional boundary.

Relevant logs and/or images

Environment information

• OpenFOAM version :v2206

Possible fixes

I thought that the bug was caused by robin boundary condition.

Three parameters(valueFraction, refValue, refGrad) related robin boundary condition is calculated in turbulentTemperatureRadCoupledMixedFvPatchScalarField.C. After that, since these are based on temperature, it is converged to an enthalpy based value at mixedEnergyFvPatchScalarField.C.

The conversion of valueFraction seems to assume constant specific heat. When I modified like below, the heat flux at multiregional boundary can be balanced.

original :
turbulentTemperatureRadCoupledMixedFvPatchScalarField.C
valueFraction() = KDeltaNbr/(KDeltaNbr + alpha);

modified :
turbulentTemperatureRadCoupledMixedFvPatchScalarField.C
valueFraction() = KDeltaNbr/(KDeltaNbr + alpha) * myThermo->Cp(pp, Tp, patchi) *(TcNbr-Tc)/(myThermo->he(pcn, TcNbr, samplePatchi) - myThermo->he(pc, Tc, patchi));