[MchEngNrd]

Right now I have a simple simulation with several different solid regions and a 'dummy' fluid region (which has not be developed yet). Right now, I am trying to establish the best way of recording/obtaining the heat flux between regions (namely solid-to-solid at the moment).

At all of my interfaces I use turbulentTemperatureCoupledBaffleMixed since it allows for interfacial resistance to be coded in (something I need). My code runs fine, but I would like to be able to record the heat flux across the patches between regions.

The most straightforward solution is using the wallHeatFlux utility, but the calculation is simply done using:

        surfaceScalarField heatFlux
        (
            fvc::interpolate
            (
                (
                    turbulence.valid()
                  ? turbulence->alphaEff()()
                  : thermo->alpha()
                )
            )*fvc::snGrad(h)
        );

which clearly a very different calculation than the one used in turbulentTemperatureCoupledBaffleMixed.

My question is, what is the best way to obtain the values calculated by turbulentTemperatureCoupledBaffleMixed? I noticed it creates the fields:

makePatchTypeField
(
    fvPatchScalarField,
    turbulentTemperatureCoupledBaffleMixedFvPatchScalarField
);

and:

tmp<scalarField> myKDelta = kappa(*this)*patch().deltaCoeffs();

Can I somehow access one of these fields during runtime using SWAK or something similar? Or would it be easier to replicate the code used in turbulentTemperatureCoupledBaffleMixed in a new post-processing utility (i.e. "my_wallHeatFlux")? The only other option I see is creating my own "my_compressible" turbulence sub-model, which has a different version of turbulentTemperatureCoupledBaffleMixed code that includes a write function.

I'm in a little over my head and its hard to tell what fields are accessible and how to access them. If someone could point me the right direction, I'd be most appreciative. I haven't found the answer despite many different searches on this site.

[hcl734]

Did you find an answer to your question?
If so would you be so kind to share it with me

[MchEngNrd]

Not yet, but its primarily because I haven't been working on this side of my project lately. I'll update this thread if I establish anything.

[Bloerb]

I think I am confused. Please check if I am correct. Let me begin by quoting Hrvoje Jasak from here.

1) a mixed boundary condition is a combination of a fixedValue and a fixedGradient boundary condition, controlled by the valueFraction variable.

For valueFraction = 1, the mixed b.c. behaves as fixedValue; for valueFraction = 0 you get the fixedGradient. Values in between are a blend of the two.

2) The mixed b.c. is therefore specified by three fields:

- a refValue field, giving the value for the fixedValue part
- a refGrad field, giving the gradient for the fixedGradient part
- a valueFraction, telling the mixed b.c. what to do.

The turbulentTemperatureCoupledBaffleMixed boundary condition calculates the temperature by harmonic averaging.

The first step is looking up the T field of the other region: nbrField... = ...TnbrName_
Afterwards we set this into the variable nbrIntFld and calculate the value the heat conduction kappa has in the other region and devide it by the cell-to-cell distance which is the deltaCoeffs part.

    // Swap to obtain full local values of neighbour internal field
    tmp<scalarField> nbrIntFld(new scalarField(nbrField.size(), 0.0));
    tmp<scalarField> nbrKDelta(new scalarField(nbrField.size(), 0.0));

    if (contactRes_ == 0.0)
    {
        nbrIntFld() = nbrField.patchInternalField();
        nbrKDelta() = nbrField.kappa(nbrField)*nbrPatch.deltaCoeffs();
    }
    else
    {
        nbrIntFld() = nbrField;
        nbrKDelta() = contactRes_;
    }

Now we do so for the region the patch is in

 tmp<scalarField> myKDelta = kappa(*this)*patch().deltaCoeffs();

And now we set the new values of our patch.

    this->refValue() = nbrIntFld();

    this->refGrad() = 0.0;

    this->valueFraction() = nbrKDelta()/(nbrKDelta() + myKDelta());

Which is setting the ref temperature to the one of the other patch and sets the ref gradient to zero. In addition the valueFraction we get from the kappa values tells us how to average those.

A few lines later in the debug options it reads

 scalar Q = gSum(kappa(*this)*patch().magSf()*snGrad());

which is the heat transfer Q and this is essentially the same as the one calculated in the heatflux utility. If someone could give me a more detailed explanation of this boundary condition with some explanations of the code I'd be delighted, because I'm also not all that sure about everything.

You should however be able to extract the values with a codedfunctionobject ().

One example of this here:

functions
(
    pAverage
    {
        functionObjectLibs ("libutilityFunctionObjects.so");
        type coded;
        redirectType average;
        outputControl outputTime;
        code
        #{
            const volScalarField& p =
                mesh().lookupObject<volScalarField>("p");
            Info<<"p avg:" << average(p) << endl;
        #};
    }
);

[zfaraday]

Hi Stephan,

You seem to understand quite good the code behind turbulentTemperatureCoupledBaffleMixed BC, at least what you describe is ok so I don't know exactly what are your hesitations about it. Maybe, you are not totally aware of what refValue, rafGrad and valueFraction mean and how are they used in the calculations. If this is what you don't know, keep in mind that a boundary value is computed in OpenFOAM according to the expression below:



Where f is the valueFraction, Tp the value of temperature at the patch and us the inverse value of the deltaCoeffs() function. Ti is the value at the cell center.

To get a better understanding about how bounday conditions are specified in OpenFOAM you can take a look at this document, I wrote during my Master Thesis development. After reading it maybe you are able to understand better each variable defined in turbulentTemperatureCoupledBaffleMixed BC.

Hope it helps.

Besst regards,

Alex

[matteoL]

Dear all,
what I don't get is why the reference gradient is set to zero. I would have set it to

nbrField.kappa(nbrField)*nbrField.snGrad() to match the neighbor heat flux

Hope someone can help.

Matteo

[zfaraday]

Hi Matteo,

Maybe the explanation comented out in the source code itself can help you. Below you can see the whole comment:

    // Both sides agree on
    // - temperature : (myKDelta*fld + nbrKDelta*nbrFld)/(myKDelta+nbrKDelta)
    // - gradient    : (temperature-fld)*delta
    // We've got a degree of freedom in how to implement this in a mixed bc.
    // (what gradient, what fixedValue and mixing coefficient)
    // Two reasonable choices:
    // 1. specify above temperature on one side (preferentially the high side)
    //    and above gradient on the other. So this will switch between pure
    //    fixedvalue and pure fixedgradient
    // 2. specify gradient and temperature such that the equations are the
    //    same on both sides. This leads to the choice of
    //    - refGradient = zero gradient
    //    - refValue = neighbour value
    //    - mixFraction = nbrKDelta / (nbrKDelta + myKDelta())

Besides that, in this document, that I wrote during the developement of my Master's Thesis, I explain the mathematical deduction of each term (refValue, refGrad and f) of the equation that OpenFOAM uses in order to compute the boundary values. I hope that after reading my document you can figure everything out and understand why refGradient is set to 0.

Hope it helps.

Best regards,

Alex

[matteoL]

Thank you very much for your quick answer!
the document is indeed very useful.

I was trying to understand a bit more into deep the way this b.c. works because in my industrial multi solid case I have one of the solids whose temperature converges very slowly (approx. 5000 iterations, while for the other solids it takes around 1000). Have you experienced anything like this (slow solid temperature convergence)? Is there any trick to speed up solid convergence? I'm running steady simulation.
In my case solid heat balance doesn't show any imbalance (residuals go down fairly well). It's just the coupling with the fluid around it that takes ages to converge...

Best regards

[zfaraday]

The matter of the slow convergence of the chtMR solvers is a well known issue. Take a look at the thread below where some ways to speed up convergence are discussed.

conjugate-heat-transfer-slow-solid-temperature-convergence

Btw, as for the issue you mention of only one solid converging slowly, the convergence in solid regions depends on its physical properties too. Is there a big difference on its thermal conduction or density compared to the other solids?

Best regards,

Alex

[matteoL]

Thanks,
Unfortunately i'm already using UR=1 for the solids and 0.95 for the fluid.
Also for each iterations solid residuals go down very well so using more non orthogonal correction wouldn't make any difference.
Why would you expect solid density to make some difference in a steady state simulation?

Regards

[zfaraday]

What about the solver tolerances (relTol and tolerance)? Maybe you can lower them after the first iterations so as to get a better convergence (remember that in these kind of cases they should have very low values, somewhere in the forum I read a while ago that the recomended values for tolerance were around 1e-10 or even lower!).

I think you are right, I checked the code and density seems not to have any influence on convergence for steady state cases. However, it is important to notice that the values of thermal conductivity (and probably heat capacity too) have a strong influence over the convergence time in this solver. Maybe you could even give a lower value to conductivity for the first time steps and then increase it to its desired value (I don't remember if the relation between conductivity and convergence velocity goes in this direction or in the opposite direction) in order to increase convergence velocity.

Best regards

[matteoL]

Dear Alex,
can you confirm that

refGradient = zero gradient
refValue = neighbour value
mixFraction = nbrKDelta / (nbrKDelta + myKDelta())

can be derived just in the case of conformal mesh between the two sides?
in this case the areas would be the same and so cancel on both sides of heat fluxes equation. In my opinion this relationship is no longer valid in case of non conformal interfaces. What do you think?

Regards

[Bloerb]

Not quite. This is still true for non conformal meshes but the method to compute the values (marked below) needs to change.

refGradient = zero gradient
refValue = neighbour value
mixFraction = nbrKDelta / (nbrKDelta + myKDelta())

If your mesh is not identical your faces do not have the same area and orientation. This can lead to problems as you pointed out. The way refValue and mixFraction is calculated is done via sampling. In your polyMesh/boundary file you can change the sample mode used. The standard option is as follows:

sampleMode      nearestPatchFace;

This sets the temperature of each face in your patch to the nearest face of the neighbor patch. Switching to something like nearestPatchFaceAMI can increase accuracy in some situations.

IMO though you should always make sure that your meshes are identical on coupled patches. Especially coupling patches of different sizes to one another should not be done. Interpolating can drastically reduce accuracy. SnappyHexMesh as well as commercial meshing programs like ansa do have options to make sure this is the case. And should be used! A good way to estimate these loses is to calculate the wallheatflux. A big difference between the two regions can easily be just because of poor interpolation.

And your problem with slow convergence is something I had problems with when I started using these solvers. I suppose you have meshed your regions with tetrahedral elements? I experienced much quicker convergence on polyhedral meshes. If you are using ansa just use the convertToPoly button. If your mesh is already fine enough you can not see a difference between the results but convergence is progessing a lot faster from my experience. OpenFoam also offers the polyDualMesh utility to convert a tet mesh to polys.

[matteoL]

Thanks Bloerb,

I tried and having a conformal interface didn't help the convergence.
Would you transform to poly the fluid, the solid region or both?

Changing a bit topic:
If I put anisotropic conductivity in one of my solid region I get to a converged solution but with non zero heat balance!
The values of the heat fluxes at each boundary make sense (I get the same values on the solid and fluid sides). It's just the Total heat flux for the anisotropic solid region that is not zero! Something like 100W for a 1000W imposed heat rate.
Changing schemes for sngrad (to corrected) or mesh didn't help.
Does anybody have experience with anisotropic conductivities?

Thanks a lot

[Bloerb]

If i understood you correctly your heat fluxes cancel out on coupled patches as they should, but the imposed heat flux on an exterior boundary is much higher than specified?
Please post how you specified your anisotropic conductivity, schemes and solution settings. Or even better upload everything without the mesh and I'll check.

[volker1]

Hi there,

I have a comment/question related to heat conservation in the turbulentTemperatureCoupledBaffleMixed BC used by chtMultiRegionFoam.
Possibly this is related to the heat balance problem stated above.

Although this BC formally bases on equal heat flows on both sides, meanwhile I think that in practice the BC will not always operate strictly heat conservative.

Within each time step, the fluid and solid regions are not solved in a common, big matrix but sequentially region by region. In chtMRF first the loop over all fluid regions is completed, followed by the loop over all solid regions.
The first call of this BC is always done from the fluid side. Depending on kappa, mesh dimensions and time step, in the inner iterations of a time step the heat flow from the boundary may not only reach the face cell of the boundary, but also get deeper into the bulk cells of the fluid volume.
As a consequence, a certain amount of heat will already be gone into the fluid bulk when the second call for this boundary patch from the solid side is done. The heat balance of the first iteration on the solid side is equal to the heat balance done at the final iteration of the fluid side (it does not look like that the solid side uses the old fluid temperature field from the previous time step for heat balance evaluation, right?). As stated above, this means that the difference between the first and the final iteration of the previous fluid side part of the heat balance will not be taken into account.
Because a similar difference may occur on the solid side, too, some kind of compensation may take place so the actual error might be much smaller than the above amount of heat that already went into the bulk.
I think the product of kappa (thermal conductivity) and deltaCoeffs() (=inverse distance from face cell center to boundary) should be the decisive quantity for the accuracy of the heat balance by the BC. If this product is equal on both sides, heat conservation will turn out close to perfect; otherwise it might however become poor depending on the difference of the product kappa*deltaCoeffs() on the fluid and the solid side.

It would help me a lot if someone could clarify if there is an error in my considerations or if this BC in fact works as I suspect.

Regards, Volker

[Bloerb]

From my interpretation of the code this is correct.

[volker1]

Hi Bloerb,

Thanks for your confirmation.
My C++ capabilities are growing but are still limited and so I have my doubts.

My background to stumble on this issue is that I made a mass transfer BC based on the thermal BC. I used the thermal solver and BC as draft and more or less just replaced thermal conductivity by diffusion and temperature by concentration. Formally physics are quite the same and everything seemed more or less straightforward. As things started to look reasonable, I tried a quantitative mass balance. The quantitative results were less pleasant than the pretty ParaView colour graphs suggested and I tracked my problems through to the above issue.
As my solid and fluid diffusivities differ by orders of magnitude this inaccuracy probably came out more obvious than in the original thermal BC. For now it looks like there is no further need to check my mass transfer BC for another programming mistake.

Instead my present approach is to force the mass balance by calculating the mass transfer from fluid to solid only once per patch field in the fluid loop before the matrix solution and inner iterations start. The mass transfer from/to the solid neighbour is collected and written into an additional source term field that remains constant during a time step. For the corresponding call from the solid side in the same time step, the same but negative mass transfer from the fluid patch field is written into the opposite face cells of the corresponding source term field for the solid. This way a new mass transfer calculation on the solid side is avoided.
The BC on both sides just sets this->valueFraction() = 0 and this->refGrad() = 0, so it practically acts as a "zeroGradient" BC. Consequently everything related to surface mass transfer between the regions is put into the equations via the opposite sign source term fields and therefore the mass balance should come out fine.
Well, I think I will find out if it really does.

Thanks again!

Regards, Volker

[Nomitude]

Hello Matteo,

did you by any chance find out more about the strange behaviour of the tTCBM boundary when one or both sides are anisotropic?

Originally Posted by matteoL

Changing a bit topic:
If I put anisotropic conductivity in one of my solid region I get to a converged solution but with non zero heat balance!
The values of the heat fluxes at each boundary make sense (I get the same values on the solid and fluid sides). It's just the Total heat flux for the anisotropic solid region that is not zero! Something like 100W for a 1000W imposed heat rate.
Changing schemes for sngrad (to corrected) or mesh didn't help.
Does anybody have experience with anisotropic conductivities?

Thanks a lot

I'm in a bit of a struggle because I need to simulate a highly anisotropic solid with different solids (and fluids with frozenFlow) around. Everything works fine as long as all regions are isotropic, but if one region has anisotropic thermal conductivity, the heat fluxes don't add up a bit, even if I set the conductivity to the same value in all directions and all regions have effectively the same thermophysical properties. This is the thermophysicalProperties file for the region with "anisotropic" thermal conductivity:

thermoType
{
    type            heSolidThermo;
    mixture         pureMixture;
    transport       constAnIso;
    thermo          hConst;
    equationOfState rhoConst;
    specie          specie;
    energy          sensibleEnthalpy;
}

coordinateSystem{type cartesian;
origin (0 0 0);                coordinateRotation
                {
                    type    STARCDRotation;
                    rotation      (0 0 0);
                }}

mixture
{
    specie
    {
        nMoles      1;
        molWeight   12;
    }

    transport
    {
        kappa   (10 10 10);
    }

    thermodynamics
    {
        Hf      0;
        Cp      1000;
    }

    equationOfState
    {
        rho     3000;
    }
}

I would be reeeeeally grateful for every hint anybody can give me!!

Thanks,
Stephan