[hcl734]

Hi all,

I simulate a heat exchanger using chtMultiRegionSimpleFoam and OF 2.4.0.
The results are somewhat realistic but I think there might be a problem with thermal convergence in the SOLID because the heat balance doesn't sum up.

heat balance solid:

The heat transferred at water side (blue) fluctuates and tops the heat transferred from exhaust gas side (red) which makes no sense for me



Other parameters as global temperature increase, pressure drop and the residuals point to a converging solution (become constant) although the maximum temperature in the fluid rises continously.

I measure the transferred heat with a functionObject

    energyAbsorbedSOLID
    {
        type patchExpression;
        outputControlMode timeStep;
        outputInterval 1;
        accumulations (
            sum
        );
        patches (".*");
        region SOLID;
        variables ("kvar=140;"
                   );
    expression "kvar*(snGrad(T))*area()";
        verbose true;
    }

[hcl734]

Another thing are the high time step continuity errors, which decrease but are still quite high



Could somebody point me to a good definition of those values so I can understand them better?

[hcl734]

It seems to be a matter of convergence.
If I continue the run it becomes constant with a small difference.
Increasing nNonOrthCorrectors in the Solid was also helpful, since more iterations are done for the solid then, which let the case converge faster.

[Bloerb]

Increasing nNonOrthCorrectors in the Solid was also helpful, since more iterations are done for the solid then, which let the case converge faster.

hmm that is true but there is a better way...Please post your fvSchemes and fvSolution for your solid region. Heat transfer calculations need a very low residual to converge. A few tips about this in general:

Check your residuals (v2.4.0 has a neat function object for this). If there is no change check the log file if the h equation does iterate. You might want to use minIter 1 to be safe. The residuals alone however are often misleading. Always use the wallHeatFlux utility to check if the heatflux from one region to the next is the same. Then check again in paraview. If the heat flux isn't smooth your mesh might just be to coarse (always check the y+ value. I'd recommend a value below 1). Problems arise if the mesh isn't identical on bothsides which can be significant. Using prism layers helps alot. Temperature on a few patches or min/max values should also be checked. For all of those there are function objects to use.

In your case it is probably the relaxation factor.
The relaxation factor for h (or e depending on your thermoPhysicalProperties) has a huge influence on convergence speed. I only use them to stop the simulation from crashing in the first few timesteps. I'd fully remove them if possible (in both regions). Another reason might be your mesh or discretisation.

[hcl734]

Hi Stephan,

I am using OF 3.0 now, just to avoid confusion.


fvSchemes for WATER

ddtSchemes
{
    default steadyState;
}

gradSchemes
{
    default         Gauss linear;
    grad(U)         cellLimited Gauss linear 1;
}

divSchemes
{
    default         none;
    div(phi,U)      bounded Gauss upwind;
    div(phi,h)      bounded Gauss upwind;
    div(phi,k)      bounded Gauss upwind;
    div(phi,K)      bounded Gauss upwind;
    div(phi,epsilon) bounded Gauss upwind;
    div(phi,omega) bounded Gauss upwind;
    div(phi,R)      bounded Gauss upwind;
    div(R)          Gauss linear;
    div((muEff*dev2(T(grad(U))))) Gauss linear;
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear limited 0.333;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         limited 0.333;
}

fluxRequired
{
    default         no;
    p_rgh;
}

wallDist
{
    method meshWave;
}

fvSolution for WATER

solvers
{
    rho
    {
        solver          PCG
        preconditioner  DIC;
        tolerance       1e-7;
        relTol          0;
    }

    p_rgh
    {
        solver           GAMG;
        tolerance        1e-7;
        relTol           0.01;

        smoother         DIC;

        cacheAgglomeration true;
        nCellsInCoarsestLevel 3126;
        agglomerator     faceAreaPair;
        mergeLevels      1;

        maxIter          100;
    }

    "(U|h|k|epsilon|omega)"
    {
        solver           PBiCG;
        preconditioner   DILU;
        tolerance        1e-7;
        relTol           0.1;
    }
}

SIMPLE
{
    momentumPredictor on;
    nNonOrthogonalCorrectors 2;
    pRefCell        0;
    pRefValue       100000;
    rhoMin          rhoMin [1 -3 0 0 0] 1000;
    rhoMax          rhoMax [1 -3 0 0 0] 1100;
}

relaxationFactors
{
    fields
    {
        rho             1;
        p_rgh           0.5;
    }
    equations
    {
        U               0.3;
        h               0.7;
        nuTilda         0.7;
        k               0.7;
        epsilon         0.7;
        omega           0.7;
        "ILambda.*"     0.7;
    }
}


relaxationFactors
{
    fields
    {
        rho             1;
        p_rgh           0.4;
    }
    equations
    {
        U               0.5;
        h               0.9;
        nuTilda         0.9;
        k               0.9;
        epsilon         0.9;
        omega           0.9;
        "ILambda.*"     0.9;
    }
}

fvSchemes for SOLID

ddtSchemes
{
    default     steadyState;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;
}

laplacianSchemes
{
    default         none;
    laplacian(alpha,h)  Gauss linear corrected;
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         uncorrected;
}

fluxRequired
{
    default         no;
}

wallDist
{
    method meshWave;
}

fvSolution for SOLID

solvers
{
    h
    {
        solver           PCG;
        preconditioner   DIC;
        tolerance        1e-06;
        relTol           0.1;
    }
}

SIMPLE
{
    nNonOrthogonalCorrectors 10;
}

relaxationFactors
{
    fields
    {
    }
    equations
    {
        h               0.7;
    }
}

[hcl734]

My residuals go down and become constant. For h at 1e-6 but it takes sometime.
I also monitored pressure drop, temperature increase and maximum temperature in fluid and solid which become constant as well.
Heat flux I checked with a functionObject for the SOLID region

    energyAbsorbedWATER
    {
        type patchExpression;
        outputControlMode timeStep;
        outputInterval 1;
        accumulations (
            sum
        );
        patches ("SOLID_TO_WATER");
        region WATER;
        variables ("kvar=0.5;"
                   );
    expression "(snGrad(T))*area()";
        verbose true;
    }

Meshes on SOLID and WATER side don't match and I am using mappedWall and AMI to couple them.
The average yPlus is below one for all walls.

[hcl734]

I also wanted to check the heat balance in the fluid.
And to obtain the enthalpy flow difference between IN- and OUTLET I used

   energyOutlet
   {
        type patchExpression;
        outputControlMode timeStep;
        outputInterval 1;
        accumulations (
            sum
        );
        patches ("OUTLET");
        region WATER;
    variables (
                   );
    expression "(A+B*T-C*sqr(T))*T*phi";
        verbose true;
   }

witch A, B, C being the coefficients of my cp-polynomial.
But there seems to be some difference between what is going into the water and what I get from subtracting enthalpy flow for inlet and outlet.

[Bloerb]

If your residuals become constant (in one region) make sure that this is not due to not solving anymore.

Solving for h, Initial residual = .., Final residual = .., No Iterations 0

Lowering your tolerance and maybe adding a minimum number of iteration avoids this. The solver stops if the resdiual 1e-06 is reached but this is not always already fully converged. The relaxation factor has a huge impact on performance and might just block changes. Therefore I'd set them to 1.0 in all Regions. This will also speed up your simulation quite a bit. (maybe 1000 iterations instead of 10000) If it crashes you might want to increase them to 1.0 a few iterations into solving. This should be done in all regions.

solvers
{
    h
    {
        solver           PCG;
        preconditioner   DIC;
        tolerance        1e-08;
        relTol           0.1;
        minIter     1;
    }
}

SIMPLE
{
    nNonOrthogonalCorrectors 3;
}

relaxationFactors
{
    fields
    {
    }
    equations
    {
        h               1.0;
    }
}

The fvSchemes and fvSolution files however do look decent. Your y+ value is also absolutely fine. If these changes did not help I suspect some problem with the mesh. Since your mesh differs between regions interpolation losses have to be expected. One way I check if this is the case is looking if the heatlfux calculated with the wallHeatFlux utility looks decent. Often you can see dots or stripes if there is a big difference. This should look smooth like a temperature field.

Do the same problems occur whitout using polynomial transport?

[hcl734]

1) Number of iterations

Originally Posted by Bloerb

If your residuals become constant (in one region) make sure that this is not due to not solving anymore.

Solving for h, Initial residual = .., Final residual = .., No Iterations 0

Lowering your tolerance and maybe adding a minimum number of iteration avoids this. The solver stops if the resdiual 1e-06 is reached but this is not always already fully converged.

You are right, thats what happens after some time.

Solving for solid region SOLID
DICPCG:  Solving for h, Initial residual = 9.999899e-07, Final residual = 9.999899e-07, No Iterations 0
DICPCG:  Solving for h, Initial residual = 9.999774e-07, Final residual = 9.999774e-07, No Iterations 0
DICPCG:  Solving for h, Initial residual = 9.999774e-07, Final residual = 9.999774e-07, No Iterations 0
DICPCG:  Solving for h, Initial residual = 9.999774e-07, Final residual = 9.999774e-07, No Iterations 0
DICPCG:  Solving for h, Initial residual = 9.999774e-07, Final residual = 9.999774e-07, No Iterations 0
DICPCG:  Solving for h, Initial residual = 9.999774e-07, Final residual = 9.999774e-07, No Iterations 0
DICPCG:  Solving for h, Initial residual = 9.999774e-07, Final residual = 9.999774e-07, No Iterations 0
[...]

Do you know how to reduce the tolerance of residual because something like

    residualControl
    {
        p_rgh          1e-6;
        U              1e-8; 
        T              1e-6; 

        // possibly check turbulence fields
        "(k|epsilon|omega)"    1e-7; 
    }

Doesn't work according to http://www.cfd-online.com/Forums/ope...implefoam.html

Or how I can add a minimum number of iterations?

[Bloerb]

This is done as mentioned above. You just need to lower your tolerance or add minIter (again highlighted below) in your fvSolution files. This is one of the reasons heat transfer simulations should not be judged by residual alone. Much more import are the temperature values or the heat fluxes. Residual Control only aborts a calculation once the inital residual reaches a certain value. The solver on the other side stops solving once the tolerance or the relative tolerance from one time step to the next is reached . It is however correct that residual control is not implemented in this solver, maybe in part for this reason.

solvers
{
    h
    {
        ....
        tolerance        1e-08;
        relTol           0.1;
         minIter     1;
    }
}

relaxationFactors
{
    fields
    {
    }
    equations
    {
        h               1.0;
    }
}

[hcl734]

Hi,

your suggestions have been pure gold as they worked out for my case pretty good.
But now I am having another question:
It seems like I am using a significant amount of energy (17 W) when coupling solid and fluid region with compressible::turbulentTemperatureCoupledBaffleMix ed.
Some loss makes sense since I have Arbitrary Mesh Interfaces (AMI) which I couple using nearestPatchFaceAMI. But 17 W is too much and not confirmed by CFX calculations.
So I wonder whether it is a matter of thermophysicalProperties
I am using temperature dependent values for heat capacity (hPolynomial) and coupling is specified as

    {
        type            compressible::turbulentTemperatureCoupledBaffleMixed;
        value           uniform 348.15;
        Tnbr            T;
        kappa           solidThermo; //fluidThermo
        kappaName       none;
    }

Is this correct?
Is solidThermo temperature dependent or should I use the lookup entry?


Description of coupling condition

Description                                    Mixed boundary condition for temperature, to be used for heat-transfer                                    on back-to-back baffles.                                
                                   Specifies gradient and temperature such that the equations are the same                                    on both sides:                                    - refGradient = zero gradient                                    - refValue = neighbour value                                    - mixFraction = nbrKDelta / (nbrKDelta + myKDelta())                                
                                   where KDelta is heat-transfer coefficient K * deltaCoeffs                                
                                   Example usage:                                        myInterfacePatchName                                        {                                            type        compressible::turbulentTemperatureCoupledBaffleMixed;                                            neighbourFieldName  T;                                            kappa               lookup;                                            kappaName           kappa;                                            value               uniform 300;                                        }                                
                                   Needs to be on underlying mapped(Wall)FvPatch.                                
                                   Note: kappa : heat conduction at patch. Gets supplied how to lookup                                        calculate kappa:                                    - 'lookup' : lookup volScalarField (or volSymmTensorField) with name                                    - 'fluidThermo' : use fluidThermo and compressible::RASmodel to calculate                                        kappa                                    - 'solidThermo' : use solidThermo kappa()                                    - 'directionalSolidThermo' directionalKappa()                                
                                   Note: runs in parallel with arbitrary decomposition. Uses mapped                                    functionality to calculate exchange.

[Bloerb]

solidThermo is correct for solid regions and fluidThermo correct for fluid regions. This means that kappa is determined via the thermophysical model used. (In your case the polynomial ones) Since the "turbulentTemperatureCoupledBaffleMixed" boundary condition can be used with other solvers and solvers written by yourself it has the lookup functionality. lookup sets kappa to the value specified in the kappa file. (You need to create one and probably modify the chtmultiregion solvers to use it). The name of the file needs to be set by the kappaName entry. Which is why it right now is set to none. If no file is present lookup should however be the same as solidThermo.

One way to check where your looses occur is to use the mapFields utility to map one boundary field to the neighbour region and afterwards subtract them with foamCalc. With that you should see problematic areas. The syntax for this is a bit complicated but you could give it a try.

Are the results using the nearestPatchFaceAMI better than with the other options available or why did you switch? And can you post a closeup of the boundary between the regions so I can get a feeling for the differences between the mesh points of the two regions?

I do not expect your polynomial to influence your losses.