[zfaraday]

As I stated in a previous post (wrong calculation of wallheatflux utility in solid patches) something in wallHeatFlux utility does not do what it should. However, I just found something out even more disturbing...

Well, now I'm doing some simple tests with a couple of solids (actually it could have been with only one region since both solids have the same properties) with chtMultiRegionFoam. It's a cube (made of two halves of steel regions) isolated that receives a convective heat flux from one of the faces (the one on the left in the attached picture) and realeses heat from another face also by convection (the one on the right).



As you can see, a very simple case. I'm testing the validity of some custom BC's I created with groovyBC, but, first of all, I wanted to check some values with basic convection BC.

I run twice the same case with the same BC's for heat transfer in both caces, one of them using externalWallHeatFluxTemperature BC and the other one with groovyBC. As expected, both cases gave me exactly the same temperature distribution. So far, Good.

However, as I wanted to check heat balance, I added this piece of code to controlDict:

functions
{
    energyVariation  // thermal power variation from the previous timeStep to the current one
    {
        type swakExpression;
        region base1;
        outputControlMode timeStep;
        outputInterval 1;
        aliases {
            rho thermo:rho;
        }
        variables (
        "Eb2{/base2}=sum(rho*vol()*(h-oldTime(h))/deltaT());"
        );
        expression "sum(rho*vol()*(h-oldTime(h))/deltaT())+Eb2";
        accumulations ( average );
        valueType internalField;
        verbose true;
    }
    energyAbsorbedBase1  //Thermal power abosorved by all patches of base1 region
    {
        type patchExpression;
        outputControlMode timeStep;
        outputInterval 1;
        region base1;
        patches (".*");
        variables ("k=50;"
                   );
	expression "k*(snGrad(T))*area()";
        accumulations (
            sum
        );
        verbose true;
    }
    energyAbsorbedBase2  //Thermal power abosorved by all patches of base2 region
    {
        type patchExpression;
        outputControlMode timeStep;
        outputInterval 1;
        region base2;
        patches (".*");
        variables ("k=50;"
                   );
	expression "k*(snGrad(T))*area()";
        accumulations (
            sum
        );
        verbose true;
    }
    totalEnergyAbsorbed   //Total thermal power abosorved by all patches (sum)
    {
        type swakExpression;
        region base1;
        valueType patch;
        patchName minX;
        outputControlMode timeStep;
        outputInterval 1;
        aliases {rho thermo:rho;}
        variables ("k=50;"
                   "minX=sum(50*(snGrad(T))*area());"
                   "maxX{maxX/base2}=sum(50*(snGrad(T))*area());"
                   "maxZ1{maxZ}=sum(50*(snGrad(T))*area());"
                   "maxZ2{maxZ/base2}=sum(50*(snGrad(T))*area());"
                   "minZ1{minZ}=sum(50*(snGrad(T))*area());"
                   "minZ2{minZ/base2}=sum(50*(snGrad(T))*area());"
                   "maxY1{minZ}=sum(50*(snGrad(T))*area());"
                   "maxY2{minZ/base2}=sum(50*(snGrad(T))*area());"
                   "minY1{minY}=sum(50*(snGrad(T))*area());"
                   "minY2{minY/base2}=sum(50*(snGrad(T))*area());"
                   );
        expression "minX+maxX+minY1+minY2+minZ1+minZ2+maxY1+maxY2+maxZ1+maxZ2";
        accumulations ( average );
        verbose true;
    }
}

And here comes the disaster...

Let's take a look, for instance, at the results coming from time 600!

·Extract from the log file (Exactly the same results obtained in both cases)

Time = 600


Solving for solid region base1
DICPCG:  Solving for h, Initial residual = 0.072001687, Final residual = 7.424161e-07, No Iterations 7
DICPCG:  Solving for h, Initial residual = 6.710513e-07, Final residual = 6.710513e-07, No Iterations 0
Min/max T:min(T) [0 0 0 1 0 0 0] 25.000066 max(T) [0 0 0 1 0 0 0] 25.326546

Solving for solid region base2
DICPCG:  Solving for h, Initial residual = 0.07214109, Final residual = 8.1396631e-07, No Iterations 7
DICPCG:  Solving for h, Initial residual = 7.3557421e-07, Final residual = 7.3557421e-07, No Iterations 0
Min/max T:min(T) [0 0 0 1 0 0 0] 24.915762 max(T) [0 0 0 1 0 0 0] 25.000121
ExecutionTime = 0.25 s  ClockTime = 0 s

Expression energyVariation :  average=128.47154
Expression energyAbsorbedBase1 on base1_to_base2:  sum=-0.21461241
Expression energyAbsorbedBase1 on minY:  sum=0
Expression energyAbsorbedBase1 on maxY:  sum=0
Expression energyAbsorbedBase1 on maxZ:  sum=0
Expression energyAbsorbedBase1 on minX:  sum=173.21344
Expression energyAbsorbedBase1 on minZ:  sum=0
Expression energyAbsorbedBase2 on maxX:  sum=-44.723404
Expression energyAbsorbedBase2 on minY:  sum=0
Expression energyAbsorbedBase2 on base2_to_base1:  sum=0.19619027
Expression energyAbsorbedBase2 on maxY:  sum=0
Expression energyAbsorbedBase2 on maxZ:  sum=0
Expression energyAbsorbedBase2 on minZ:  sum=0
Expression totalEnergyAbsorbed :  average=128.49004  //should be equal to the value given by energyVariation
Region: base1 Diffusion Number mean: 0.027777778 max: 0.027777778
Region: base2 Diffusion Number mean: 0.027777778 max: 0.027777778

·Extract from wallHeatFlux, for base1 region and externalWallHeatFluxTemperature BC used

Time = 600
Selecting thermodynamics package 
{
    type            heSolidThermo;
    mixture         pureMixture;
    transport       constIso;
    thermo          hConst;
    equationOfState rhoConst;
    specie          specie;
    energy          sensibleEnthalpy;
}


Wall heat fluxes [W]
maxY 0
minX 111.68017
minY 0
minZ 0
maxZ 0
base1_to_base2 -0.325975

·Extract from wallHeatFlux, for base1 region and groovyBC used

Time = 600
Selecting thermodynamics package 
{
    type            heSolidThermo;
    mixture         pureMixture;
    transport       constIso;
    thermo          hConst;
    equationOfState rhoConst;
    specie          specie;
    energy          sensibleEnthalpy;
}


Wall heat fluxes [W]
maxY 0
minX 173.21344
minY 0
minZ 0
maxZ 0
base1_to_base2 -0.325975

·Extract from wallHeatFlux, for base2 region and groovyBC used

Time = 600
Selecting thermodynamics package 
{
    type            heSolidThermo;
    mixture         pureMixture;
    transport       constIso;
    thermo          hConst;
    equationOfState rhoConst;
    specie          specie;
    energy          sensibleEnthalpy;
}


Wall heat fluxes [W]
maxY 0
maxX -44.723405
minY 0
minZ 0
maxZ 0
base2_to_base1 0.176

Let's analyze the results calmly.

First of all, the most disturbing values that can be higlighted are the values of the heat flux at minX patch (the one that recevies convective heat). wallHF utility shows a value of minX 111.68017 in the case where externalWallHFTemperature BC is used, while the value shown when groovyBC is used is minX 173.21344. This last one is exactly the same value that can be seen in the output log file calculated with the function objects.

So far, it can be stated that something goes wrong. If we take the values given by the function objects as correct, then wallHF utility is not working as it should. I have to say that I got specially surprised by these results, since the last time I tried to compare somehow both methods to compute convective HF (groovyBC and externalWallHFTemperature BC) I would say (maybe not too loud, I should have to check it again) that wallHF utility and function objects gave the same results for the convective HF at walls... The only incorrect value I found was with the coupling of regions as I stated in the post mentioned above.

Here, in the presented case, it can also be noticed the mismatch between the values given by wallHF utility and the function objects in the case of coupled patches. You can see that the output log file shows a value of base1_to_base2: sum=-0.21461241, while wallHF utility shows base1_to_base2 -0.325975 instead. Quite an important difference that I don't understand where it comes from.

Another point that may deserve to be pointed out is the difference between the outgoing HF of the region base1 and the incoming HF of the region base2 (that is to say, the HF through the coupled patch). According to the laws of thermodynamics it should be the same, otherwise some amount of energy is destroyed or created in the interregion. This is what it is happening here, nor with wallHF utility or with the function objects the balance is achieved. The difference is bigger with wallHF utility and keeps growing timeStep after timeStep. It also grows with function objects but slowlier. This difference makes the energy balance unbalanced in the overall system, energyVariation != totalEnergyAbsorbed, (see function specification in my controlDict) unless in the totalEnergyAbsorbed function object internal HF's are taken into account (it can be seen that in the specification I didn't take them into account, since it make no sense physically). However, if internal HF's were taken into account then the equation energyVariation = totalEnergyAbsorbed would be accomplished. Thus, something is not working properly in the interregion, and it seems that wallHF utility is not very capable of giving feasible results...

I will keep testing to see if I can figure something out. In the meanwhile it would be great if someone could give a quick glance at the way wallHF works to see if it is its fault or swak4foam's fault instead... well, it can also be my fault if I did something wrong...

If you have really read it all so far without falling asleep, you deserve a lollipop!

Many thanks in advance!

Regards,

Alex


Note: HF has been used as Heat Flux to shorten it and make it (a little) faster to wright.