chtMultiRegionSimpleFoam - of2112 does not catch analitical results where of2106 does

Summary

The same test case if it is run with of2106 the analytical solution is obtained, if it is run with of2112 the analytical solution is not catched.

Steps to reproduce

The case is shared in this repo. The test case is taken from this article. Just run "Allrun-serial" to launch the simulation. Try to run the simulation loading of2106 environment and then using of2112 environment. The solution given by of2016 is correct the one computed using of2112 is wrong by some degrees in temperature.

Example case

The test case is taken from this article.

h is the convective coefficient on the far right of the system, it is so low as it must represent vacuum. Lc and kc are set to 1e-6 m and 1e-4 W/m/K.

The case is very simple and fast.

What is the current bug behaviour?

Analitical solution not catched by of2112 but catched by of2106. The one computed by of2112 is 330 K and 335 K whilst the correct one is 327 K 338 K.

What is the expected correct behavior?

To catch the analytical solution.

Relevant logs and/or images

Please see section above.

Environment information

• OpenFOAM version : v2112 | v2106
• Operating system : ubuntu
• Hardware info :
• Compiler : precompiled versions of v2112 and v2106

Possible fixes

added community label

@Mattijs @Sergio

Hi, Thanks. Any idea where the problem might be: Radiation, BC's, solver?

Hi @Sergio, I think is in the BC compressible::turbulentTemperatureRadCoupledMixed because in of2106 using compressible::turbulentTemperatureCoupledBaffleMixed provided the correct analitical solution. The latter BC has been superseded by compressible::turbulentTemperatureRadCoupledMixed in of2112.

Furthermore, in the case I shared, I used compressible::turbulentTemperatureRadCoupledMixed and if I run using of2106 the analytical solution is caught but when using of2112 it does not compute the exact analytical solution. Thus the same BC works in of2106 but does not work in of2112, I think the problem is somewhere there.

Each version has its own solver and BC. Can you run the case using v2112 solver but with the BC from of2106 ?. You'll need to rename the BC into something like turbulentTemperatureRadCoupledMixed_2106

Can you run the case using v2112 solver but with the BC from of2106 ?

I was updating the comment above about this. Yes I can but it provides the same wrong results as if I run with compressible::turbulentTemperatureRadCoupledMixed.

I can still use compressible::turbulentTemperatureCoupledBaffleMixed with of2112, it just appear the following warning:

--> FOAM Warning : 
    From Foam::compressible::turbulentTemperatureCoupledBaffleMixedFvPatchScalarField::turbulentTemperatureCoupledBaffleMixedFvPatchScalarField(const Foam::fvPatch&, const Foam::DimensionedField<double, Foam::volMesh>&, const Foam::dictionary&)
    in file turbulentFluidThermoModels/derivedFvPatchFields/turbulentTemperatureCoupledBaffleMixed/turbulentTemperatureCoupledBaffleMixedFvPatchScalarField.C at line 126
    This BC has been superseded by compressible::turbulentTemperatureRadCoupledMixed which has more functionalities and it can handle the assemble coupled option for energy. 
    Adding to radiations

Are you using a kappaLayer and thickness Layer on the turbulentTemperatureRadCoupledMixed BC ?.

Yes I do use kappaLayer and thicknessLayer

Compare both cases WITHOUT kappaLayer and check if results from different versions are the same

I've just test it, they do provide the same results, does that mean that in of2112 'compressible::turbulentTemperatureRadCoupledMixed' is not using/supporting the use of kappaLayer?

It seems that in of2112 the kappalayer info are not considered while in of2106 they are taken into account.

https://www.openfoam.com/news/main-news/openfoam-v2112/boundary-conditions.

Deeply sorry, I was jumping between of2106 and of2112 to try to solve this other issue #2387. Didn't check the changelog, my bad.

In any case, I do not understand why the following syntax as reported in the documentation is not working

    {
       type            compressible::turbulentTemperatureRadCoupledMixed;
       Tnbr            T;
       kappaMethod     solidThermo;
       thicknessLayers (1e-6);
       kappaLayers     (1e-4);
       value           uniform 300;
    }

The simulation is running but the kappaLayers and thicknessLayers information is not taken into account.

@Sergio just an update: using this $FOAM_TUTORIALS/heatTransfer/chtMultiRegionFoam/multiRegionHeater with of2112, if kappaLayer and thickenssLayers are removed from the subdict in the changeDictionary file in the leftSolid and heater region the solution computed is identical as in the unmodified case.

assigned to @Sergio

I do see differences (only ran for one time step though but had to set the write precision to 16 digits). The major difference between the old and new is that the kappaLayers/Thickness can be specified for a single side separately. On this case both sides have additional resistance so if the heater-to-leftSolid is the dominant driver it will just delay the reaching of steady-state. We probably need to make the resistance one-sided for this case.

@Mattijs , note that this bug report is very closely related to #2344 (closed) The conclusion of that issue is also "The simulation is running but the kappaLayers and thicknessLayers information is not taken into account".

Based on your explaination I conclude that you would expect the resistance to be accounted for 2x, If that's the case, I really wonder why this issue and 2344 are marked as 'community' and not as 'bug' because there is clearly something wrong.

@Mattijs I am also observing odd behaviors using turbulentTemperatureRadCoupledMixed when kappaLayer is provided on both sides in 2112. I would suggest you guys to take a closer look based on the provided example above.

There is a problem due to the lack of inertia in the added resistance. Currently the boundary value is usually dominated by the non-resistance side but then only the value enters the equation. We need to re-think this - maybe enter the added resistance not only in the value calculation but in the thermophysical properties (i.e. a patch-based modification of the properties) such that the whole simulation 'sees' it.

@Mattijs Is it an option to revert the turbulentTemperatureRadCoupledMixed boundary condition in v2112 back to version v2012? Or will this give conflicts due to other changes in the code? Because now this issue is preventing us from upgrading to v2112 and I suspect it is the same for others as well.

added bug label

The v2012 implementation had its own problems - it was dependent on the old-time face values so might never converge - when setting thickness to infinitely small it still converged to a different solution compared to having no thickness.

Under investigation.

Please, test this updated BC.:

turbulentTemperatureRadCoupledMixedFvPatchScalarField.C

It seems to work for simple test cases on my side. The contact resistance needs to be specified on both sides

Thanks

mentioned in issue #2344 (closed)

@Sergio: For the updated BC using the "turbulentTemperatureRadCoupledMixedFvPatchScalarField.C" posted above, Does the contact resistance needs to be the same on both sides or can they be different? Also starting from line 400, "thicknessLayer" is not used for kDeltaC. Are you trying to switch back to 2012 or ...?

 if (thicknessLayer_ || thicknessLayers_.size())
{
    // Harmonic averaging
    {
        KDeltaC = 0.0;

        if (thicknessLayers_.size())
        {
            forAll(thicknessLayers_, iLayer)
            {
                KDeltaC += thicknessLayers_[iLayer]/kappaLayers_[iLayer];
            }
        }
        if (thicknessLayers_.size())
        {
            forAll(thicknessLayers_, iLayer)
            {
                KDeltaC += thicknessLayers_[iLayer]/kappaLayers_[iLayer];
            }
        }
        KDeltaC = 1.0/(KDeltaC + SMALL);
    }
}

@Sergio Thanks for the updated code. As @agegnehu1 pointed out in his post, I think there is a small copy-past error in the code snippet above (lines 401-422). With the correction below, the code reproduces our analytical test cases:

    if (thicknessLayer_ || thicknessLayers_.size())
    {
        // Harmonic averaging
        {
            KDeltaC = 0.0;
            if (thicknessLayer_)
            {
                const scalar t = db().time().timeOutputValue();
                 KDeltaC +=
                     thicknessLayer_().value(t)
                    /kappaLayer_().value(t);
            }
            if (thicknessLayers_.size())
            {
                forAll(thicknessLayers_, iLayer)
                {
                    KDeltaC += thicknessLayers_[iLayer]/kappaLayers_[iLayer];
                }
            }
            KDeltaC = 1.0/(KDeltaC + VSMALL);
        }
    }

@agegnehu1 , @lieven Thanks for testing and checking the code. I'll change the copy/paste error to work with both thicknessLayer_ and thicknessLayers_. The formulation is different from both previous versions. The resistance needs to be specified on both sides. The resistance is inserted between the two walls and it is added to the local gradient of the local patch of the mixed BC. This formulation is symmetric for both sides. It works well for very small resistance and large as well. I tested with the multi-word tutorial and it seems to work. As well, for simple 1D heat transfer between fixed T's on two solids.

@Sergio: I will check its response to some of my tests and let you know. When you say, "The resistance needs to be specified on both sides", do they need to be exactly the same or can they be different values on each side?

@agegnehu1 ,

The same resistance needs to be specified on both sides.