[Davide95]

Hey all!
I am a complete newbie and I'm having some troubles with the simulation of
a countercorrent heat exchanger (turbulent), stationary solution.

There are two coassial pipes. In the outer one flows hot C02 and in the inner one flows colder air.
My problem is that at the end of the simulation I have a wrong power balance: the air absorbs way more heat than the C02 releases.
(I mean that cp_air*ΔT_air*ṁ_air is not equal to cp_C02*ΔT_CO2*ṁ_CO2 ).
I have already checked that the geometry is fine, and that the cp of both are correct.
At first i thought that maybe the flux of air dissipated kinetic energy in thermal energy but looking at the solution given by openFoam I noticed that the velocity keeps almost constant .
For the solution I used chtMultiResgionSimpleFoam and a K-omega model.
If anyone can help I'm really freaking out cause seems that air takes that energy from nowhere
Thanks a lot in advance

[Bloerb]

Calculate the wall heat flux in all regions and report back. Depending on your OF version with the wallHeatFlux utility or the function object. Do they match?
Plot the temperature at a patch or the min max value over the iterations using a function object. Has this reached steady state? Post your inlet and outlet boundary conditions. How is your y+ value etc. Some examples below:
Depending on your version of OF you need to modify these. Look in the .H files inside the src/functionObjects or src/postProcessing/functionObjects folders for examples.

    Average_temp_outlet
    {
        type            surfaceFieldValue;
        libs            ("libfieldFunctionObjects.so");
        log             true;
        writeControl    timeStep;
        writeFields     false;
        surfaceFormat   none;
        regionType      patch;
        name            outlet;
        region          fluidA;
        operation       areaAverage;
        fields
        (
            T
        );
    }
    yPlus_fluidB
    {
        type            yPlus;
        region          fluidB;
        writeControl    outputTime;
        libs            ("libfieldFunctionObjects.so");
    }
    residualsFluidA
    {
      type            residuals;
      libs            ("libutilityFunctionObjects.so");

      writeControl    timeStep;
      writeInterval   1;
      region          fluidA;

      fields (p_rgh U h k omega);
    }

    residualsFluidB
    {
      type            residuals;
      libs            ("libutilityFunctionObjects.so");

      writeControl    timeStep;
      writeInterval   1;
      region          fluidB;

      fields (p_rgh U h k omega);
    }

    residualsSolid
    {
      type            residuals;
      libs            ("libutilityFunctionObjects.so");

      writeControl    timeStep;
      writeInterval   1;
      region          solid;

      fields (h);
    }

[Davide95]

When i said I am a newbie I really meant that
I tried to have a look on how to calculate the wallHeatFlux but can't really figure out.
About the inlet and outlet conditions which ones do you need (p_rgh, T, U, nut, alphat , all of them?) .
Regarding the steady state, how can I know if it has reached it? I gave it for granted but from your answer obviously is not.

Sorry for the stupid questions and thanks a lot for your patience

[Bloerb]

Lets start from the top.

You do need appropriate boundary conditions:

U

    inlet
    {
        type            flowRateInletVelocity;
        massFlowRate    0.100; // kg/s - volumetricFlowRate is also possible
        extrapolateProfile  no; // turn to yes for a developed velocity profile 
        value           $internalField;
    }
    outlet
    {
        type            pressureInletOutletVelocity;
        inletValue     uniform (0 0 0);
        value           $internalField;
    }
    fluidA_to_solid
    {
        type            noSlip;
        value           uniform (0 0 0);
    }

p

    inlet
    {
        type            calculated;
        value           $internalField;
    }
    outlet
    {
        type            calculated;
        value           $internalField;
    }
    fluidA_to_solid
    {
        type            calculated;
        value           $internalField;
    }

p_rgh

internalField   uniform 1e5;

boundaryField
{
    inlet
    {
        type            zeroGradient;
        value           $internalField;
    }
    outlet
    {
        type            fixedValue;
        value           $internalField;
    }
    fluidA_to_solid
    {
        type            fixedFluxPressure;
        value           $internalField;
    }
}

T

internalField   uniform 313.15;

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform 313.15;
    }
    outlet
    {
        type            inletOutlet;
        inletValue      $internalField;
        value           $internalField;
    }
    fluidA_to_solid
    {
        type            compressible::turbulentTemperatureCoupledBaffleMixed;
        value           $internalField;
        Tnbr            T;
        kappaMethod     fluidThermo;
    }
}

I won't bother with boundary conditions for k, epsilon, omega, alphat and nut for now. Since these are always tricky and can crash your simulation very fast. Hence always check everything WITHOUT turbulence in your first run. Once the results look plausible add turbulence. Here the kOmegaSST model should be your choice since it is arguably the most reliable and accurate. The guys doing external aerodynamics might proclaim otherwise, but we are talking about internal geometries and it is the industries choice! Most notibly because it won't crash as easily if your inlet turbulence guess is totally wrong.
All walls should have wall functions in those boundary conditions. All outlets are zeroGradient, inletOutlet, or calculated. The inlets fixedValue or calculated.

The next step is solver settings. I always advice not relaxing h or e in all regions. Introducing relaxation makes your simulation seem more stable, but it converges much much slower and crashes are just later in the run.

The next step is to check your mesh:

checkMesh -allTopology -allGeometry -region fluidA
checkMesh -allTopology -allGeometry -region solid
checkMesh -allTopology -allGeometry -region fluidB

If this fails you better have confidence in your mesh. If nonOrthogonality is above 45 add nonOrthogonal Correctors (1-5 should do the trick) If it is in the 60+-75 range you should consider remeshing. Now take a look at your fvSchemes files for the fluid. All div Schemes should be upwind. Fancy second order shit is something you can try once upwind successfully ran. Especially with compressible and heat transfer applications you need to make sure your values are bounded. If your lowest temperature inlet is 300K and you get 295K due to your second order scheme somewhere you probably shouldn't have used them in the first place. Even upwind can crash if your mesh is not perfect. Here limiters are a good choice. A quick example below:

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

laplacianSchemes
{
    default         Gauss linear limited 0.333;
}

snGradSchemes
{
    default         limited 0.333;
}

Add an fvOption file with a temperature limiter in each region to make sure it does not crash due to nonphysical temperatures in the first few iterations.

    limitT
    {
        type            limitTemperature;
        active          yes;

        limitTemperatureCoeffs
        {
            selectionMode   all;
            min             300;
            max             500;
        }
    }

Inside your controlDict add the functionObjects I have shown above. You should plot the residuals to monitor convergence. foamMonitor -l residuals.dat is one option. Once the residuals have dropped low enough (1e-6 for example) you should also check the output of the average temperature function objects. They should become stationary! Afterwards run

wallHeatFlux -region fluidA
wallHeatFlux -region fluidB
wallHeatFlux -region solid

to check if these are identical with opposite flux on each boundary. If your simulation has converged well enough this should be the case. All of these things might be slightly different in different OpenFOAM versions.
The tricky part about cht calculations is that the residuals are not always the full picture. (To be honest the aren't in all other calculations as well, but here people notice it when their results are still wrong). Always check temperature, velocity, pressure etc over the iterations next to the residuals to make sure they are not oscillating or steadily dropping.

Also calculate your y+ value. CHT calulations need a finely resolved boundary layer. Calculate your Prandtl, Reynolds, Biot, Peclet numbers and estimate the height of your boundary layer. If all the heat transfer happens in 0.1 mm and your first layer is 5cm big you won't get good results.

[Davide95]

Thanks a lot!
1) MESH
The problem is a project from the university and the mesh has been provided by the professor so I guess it's correct. I checked anyway and the Non-orthogonal factor is 37 and the y+ seems good enough.

2) Boundary conditions

Even if not identical the boundary are pretty similar to the ones you provided me. While doing the simulation I used a test case that we did in class for a simple pipe.

3) Heatflux

The heatfluxes are really different for the different walls but I saw that the difference is decreasing (even if not a lot) with the number of iterations.
I did 5000 iterations : should I do more? Someone suggested me 50k

4) THANKS A LOT
you have been super kind and it has already been really useful.

[Bloerb]

The non-orthogonality is excellent. I suppose it is a hex mesh. The most likely reason for the slow convergence is the mentioned relaxation.
Check your fvSolution files of each region. Increase your relaxation parameter to 1.0 for h or e. If it crashes lower it to 0.999 or start with 0.7 and increase it as the simulation goes on. With this 50k Iterations should not be needed. Probably not even one thousand. (This does depend on your domain, cell count and flow complexity though). A pipe heat exchanger however should have forced convection flow of high reynolds numbers into only one direction. Monitor the residuals of the flow with the function object or by monitoring the output. Make sure that in fvSolution a minIter flag is set inside each solver setting. This will prevent that the flow stops solving while the temperature has not yet converged.

"(U|k|...)"
{
...
minIter 1;
...
}

Since the region coupling is done explicitly your heatfluxes won't be exactly identical. Mapping errors due to different meshes in the regions etc simply add up. You should however come fairly close.

[Davide95]

YAY!!! I took the result of the 5k iteration and renamed it as 0. Then changed the relaxation parameter as you said and run it again! Now the heat fluxes are slighlty different but way better than before!! Thanks a lot mate!
One last question, can you explain me why changing the relaxation factor helped and how did you know that i had to change exaclty the 'h' relaxation factor?

[Bloerb]

You are solving three equations

1. momemum equation (U,p)
2. mass balance (U)
3. energy equation (h or e)

The unknowns are U,p and h. The momentum or navier stokes equation is actually a vector equation in x y and z direction, hence 3 equations, however 4 unknowns. u v w and p. The SIMPLE method solves this by adding the mass balance into the mix. The SIMPLE method also needs relaxation factors to converge. You can read up on the differences between SIMPLE and PISO to know why and what it is relaxation actually does and how the SIMPLE method functions. To put it in a nutshell: Relaxation factors make sure your solution does not explode by making the changes from one time step to the next smaller. This also effects errors, which is the good thing. as a result of this gain in stability it naturally converges slower. Remember we make changes from one iteration to the next smaller. The main reason this is necessary is due to the way the momentum and mass conservation are coupled in the SIMPLE algorithm and the non linear nature of the navier stokes equation (with the non linearity in the convective transport term). Over the years standard relaxation factors for u and p emerged. In cases with low convective transport for example the equations are "more linear" and you could increase those.

Now the energy equation (solved for h or e, which is transformed into the temperature) is linear. Hence it is much more stable and does not need the relaxation (at least in most cases). Hence changing the relaxation for h makes sure your get your results faster. The low numbers in the openfoam tutorials is a mystery to me.

[Davide95]

Super clear!
Thanks for all! Problem solved