CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)
-   -   Energy balance in solid using chtMRSF (https://www.cfd-online.com/Forums/openfoam-solving/165367-energy-balance-solid-using-chtmrsf.html)

hcl734 January 15, 2016 07:42
Energy balance in solid using chtMRSF
 
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
http://s24.postimg.org/606fppoqd/resq.png


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

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

hcl734 January 15, 2016 07:47
Another thing are the high time step continuity errors, which decrease but are still quite high

http://s12.postimg.org/l5kezoi71/Scr...5_13_45_37.png

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

hcl734 January 26, 2016 07:47
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 January 26, 2016 10:44
Quote:
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 January 27, 2016 02:58
Hi Stephan,

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


fvSchemes for WATER

Code:
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

Code:
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

Code:
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
Code:
solvers
{
    h
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance        1e-06;
        relTol          0.1;
    }
}

SIMPLE
{
    nNonOrthogonalCorrectors 10;
}

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

hcl734 January 27, 2016 03:11
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

Code:
    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 January 27, 2016 03:15
I also wanted to check the heat balance in the fluid.
And to obtain the enthalpy flow difference between IN- and OUTLET I used

Code:
  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 January 27, 2016 11:31
If your residuals become constant (in one region) make sure that this is not due to not solving anymore.

Code:
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.
Code:
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 January 28, 2016 03:10
1) Number of iterations

Quote:
Originally Posted by Bloerb (Post 582661)
If your residuals become constant (in one region) make sure that this is not due to not solving anymore.

Code:
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.

Code:
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

Code:
    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 January 28, 2016 12:24
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.
Code:
solvers
{
    h
    {
        ....
        tolerance        1e-08;
        relTol          0.1;
        minIter    1;
    }
}

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

hcl734 February 8, 2016 11:17
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
Code:
    {
        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

Code:
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 February 8, 2016 14:03
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.


All times are GMT -4. The time now is 19:17.