[keitaro7_14]

Hi,

I noticed a very strange behaviour in chtMultiRegionFoam. I prefer to share my question to the community before report a bug, because I could be wrong.

I have two solids (two squares) in contact. Very very simple. I want to know Temperature at t=200s.

All boundaries are zeroGradient except the right face, with fixed temperature.

Heat transfer is governed by the Heat Equation, where diffusivity is the key.

Solid 1 | Solid 2



Solid 1:

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  8
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      thermophysicalProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

thermoType
{
    type            heSolidThermo;
    mixture         pureMixture;
    transport       constIso;
    thermo          eConst;
    equationOfState rhoConst;
    specie          specie;
    energy          sensibleInternalEnergy;
}

mixture
{
    specie
    {
        nMoles      1;
        molWeight   100.02;   // [g/mol]
    }

    transport
    {
        kappa   10;        // [W/m/K]
    }

    thermodynamics
    {
        Hf      0;
        Cv      1000;        // [J/kg/K]
    }

    equationOfState
    {
        rho     2200;       // [kg/m^3]
    }
}

// ************************************************************************* //

Solid 2:

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  8
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      thermophysicalProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

thermoType
{
    type            heSolidThermo;
    mixture         pureMixture;
    transport       constIso;
    thermo          eConst;
    equationOfState rhoConst;
    specie          specie;
    energy          sensibleInternalEnergy;
}

mixture
{
    specie
    {
        nMoles      1;
        molWeight   100.02;   // [g/mol]
    }

    transport
    {
        kappa   10;        // [W/m/K]
    }

    thermodynamics
    {
        Hf      0;
        Cv      1000;        // [J/kg/K]
    }

    equationOfState
    {
        rho     2200;       // [kg/m^3]
    }
}

// ************************************************************************* //

Well, the result seems OK. I used an analytical expression to solve this problem for spherical coordinates, and it fits very well (using a mesh of two concentric spheres).

For this problem (regarding to the squares), I get this result from the log, when solver finishes.

Solving for solid region capa_2
DICPCG:  Solving for e, Initial residual = 0.0037510124, Final residual = 6.6981383e-07, No Iterations 3
Min/max T:333.83288 453
ExecutionTime = 0.27 s  ClockTime = 0 s

Now, I am going two change both, k and Cv in Solid 1, with the following values:
Solid 1:

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  8
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      thermophysicalProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

thermoType
{
    type            heSolidThermo;
    mixture         pureMixture;
    transport       constIso;
    thermo          eConst;
    equationOfState rhoConst;
    specie          specie;
    energy          sensibleInternalEnergy;
}

mixture
{
    specie
    {
        nMoles      1;
        molWeight   100.02;   // [g/mol]
    }

    transport
    {
        kappa   0.01;        // [W/m/K]
    }

    thermodynamics
    {
        Hf      0;
        Cv      1;        // [J/kg/K]
    }

    equationOfState
    {
        rho     2200;       // [kg/m^3]
    }
}

// ************************************************************************* //

Noticed that the thermal diffusivity is the same, so I should expect the same result. Now, the result of chtMultiRegionFoam when it finishes:

Solving for solid region capa_1
DICPCG:  Solving for e, Initial residual = 0.015590819, Final residual = 1.955389e-07, No Iterations 4
Min/max T:309.80765 371.48362

Solving for solid region capa_2
DICPCG:  Solving for e, Initial residual = 0.0090723273, Final residual = 8.2453758e-08, No Iterations 4
Min/max T:371.48425 453
ExecutionTime = 0.29 s  ClockTime = 0 s

The minimum temperature of the second block is changing!! Why? I tried to change de ddtScheme, size of the mesh, adjustTimeStep, the maxDi...impossible. I don't know if I am stupid or really there is a bug in OpenFOAM



I attach the case for OF-8, so you can test to change k and Cv mantaining the same thermal diffusivity.

Thank you so much!!

[keitaro7_14]

Well, as expected, I am the stupid here, and not the OpenFOAM developers. The problem here is that the problem is wrongly formulated. Let me explain.

When we have two solids in contact, we can model the problem like two solids of different properties in contact, or like a unique solid with variable thermal properties. It is fine, but we have to pay attention to the formulation.

I assumed that two solids can be modelled as a solid with variable diffusivity, and the equation governing his physical behaviours is the following:



but it is wrong! Why? Well,the success of the analysis/modelling depends on many factors, but the basis is the correct application ofthe heat conduction theory which is more important in case of transient thermal processes. I paid attention to the tool (OpenFOAM), but not to the formulation.

Most of the cases we can assume that there is no movement, no force field at the surface, no volume heat generation, and the pressure is constant. Neither the density nor the thermal capacity are changing with the time. Taking into account these conditions, the simplified heat conduction equation will be as it follows



Now, I can do a little transformation:



But..Can I introduce into the derivative? NO! Because they are not constant! My real model are two solids, with different thermal properties. If I assume one solid, , and change with the position. The correct formulation following this approach of one solid with variable thermal diffusivity should be as follows:



And well...we should derivate new ODEs to solve the problem.

The other approach, with two solids, is the approach of OpenFOAM, and if we have two solids, we have a system of equations.




and in the interface, we have





where



I hope this will help other researchers who have the same doubts as me and the same errors. I also recommend consulting the following publication: https://doi.org/10.1007/s10973-018-7014-4