[dvcauwe]

Dear Foamers,

I am currently working on evaluation of the heating characteristics and pressure drops in different pipe geometries. I would like to do this by simulating a short part of the tube and applying cyclic streamwise boundary conditions to obtain fully developed flow.

For this I am using a modified version of channelFoam in which the temperature equation is solved as well. However, I'm unsure on how to make the temperature cyclic. So far I have used a fixedGradient at the wall to impose constant heat flux, while scaling the temperature everywhere so that the mass-weighted averaged temperature equals the inlet bulk temperature. Is this a correct approach to obtain a developed temperature profile?

In literature I keep seeing people use dimensionless temperatures and energy source terms but I'm quite clueless on how to implement this in the solver... Apparently the basic article on this matter is by Patankar et al. (1977), “Fully developed flow and heat transfer in ducts having streamwise-periodic variation of cross-sectional area” but I can't seem to get my hands on that one for now.

I know the question has been asked before but it was never fully resolved so I would really appreciate any input/hints you might be able to offer!

Best regards,
David

[dvcauwe]

I noticed that the way Fluent does this is by scaling the entire "outlet" temperature profile with (Twall-Tbulk,out)/(Twall-Tbulk,in). This quantity is actually very useful for a constant wall temperature as it allows fast evaluation of the convection coefficient. Attached is a cylinder simulation of what I'm trying to achieve.

Now my question is how to implement this in openFOAM? I have tried using groovyBC to copy the scaled outlet values to the inlet but I think this messed up my gradients as the TEquation requires more and more iterations.

The boundary conditions in the 0/T file are as such:

boundaryField
{
wall
{
type fixedValue;
value uniform 350;
}
periodic_half0
{
type groovyBC;
patchType cyclic;
valueExpression "350-(350-TCyc)*0.995"; //0.5% increase of temperature predicted
variables "TCyc{periodic_half1}=T;";
}
periodic_half1
{
type cyclic;
}
}

Could anyone give me some hints on what I'm doing wrong and/or if there's an easier way to transfer conditions from one patch to another? jumpCyclic looks interesting also but so far I have not been able to find an example of how that works...

Best regards,

David

[mra-cfd]

Hi David,

Have you solved the problem of cyclic temperature with fixedGradient at walls? Is there a way to set the "upstream bulk temperature" like FLUENT?

Thanks,
Mohammadreza

[dvcauwe]

Originally Posted by mra-cfd

Hi David,

Have you solved the problem of cyclic temperature with fixedGradient at walls? Is there a way to set the "upstream bulk temperature" like FLUENT?

Thanks,
Mohammadreza

In this thread it's explained in a bit more detail. Basically you add an additional source term to your equations to compensate for the heat loss/gain that would normally occur in order to make your temperature field periodic.

In the case of a constant cross-section this is very easily implemented, otherwise you need to solve for lambda first. Make sure to check out the paper by Patankar about this.

Best regards,
David

[hrvig]

Originally Posted by dvcauwe

In this thread it's explained in a bit more detail. Basically you add an additional source term to your equations to compensate for the heat loss/gain that would normally occur in order to make your temperature field periodic.

In the case of a constant cross-section this is very easily implemented, otherwise you need to solve for lambda first. Make sure to check out the paper by Patankar about this.

Best regards,
David

Hi,

Did you manage to implement the source terms suggested by Patankar or did you scale the whole temperature field to match the desired upstream bulk temperature?

If you got it to work, I will be very interested to hear more

Best,
Jakob

[dvcauwe]

Hello Jakob,

In fact I never really managed to get the solving for lambda right, but most of my cases have a rotating but constant cross-section so I could just use a fixed heat flux BC. For a fully developed flow that means that you can calculate the linear temperature gradient gamma [K/m] for that heat flux and split the real temperature field into this linear gradient and a periodic field: T = Tper + gamma*x. Substituting this term into the temperature equation gives you an energy sink term -Ux*gamma on the right hand side. As in OpenFOAM performing periodic calculations including heat transfer will force you to write your own solver anyway, this solution is much easier than fiddling around with boundary conditions

Best regards,
David

[hrvig]

Originally Posted by dvcauwe

Hello Jakob,

In fact I never really managed to get the solving for lambda right, but most of my cases have a rotating but constant cross-section so I could just use a fixed heat flux BC. For a fully developed flow that means that you can calculate the linear temperature gradient gamma [K/m] for that heat flux and split the real temperature field into this linear gradient and a periodic field: T = Tper + gamma*x. Substituting this term into the temperature equation gives you an energy sink term -Ux*gamma on the right hand side. As in OpenFOAM performing periodic calculations including heat transfer will force you to write your own solver anyway, this solution is much easier than fiddling around with boundary conditions

Best regards,
David

Thanks for getting back with useful information David!

Did you implement this using the cyclic boundary conditions on inlet and outlet, just like using channelFoam in old OpenFOAM versions?

So far I am using a pimpleFoam where temperature is treated as a passive scalar. My geometry is a simple twisted pipe with constant diameter.

Best,
Jakob

[hrvig]

Hi again,

Now I implemented the source terms you suggested. I begin to believe that this method will work and I now have the cyclic temperature profile - hurra :-)

Only problem left is that I loss energy. How do you set boundaries?
Here is an overview of my solver for temperature:

fvScalarMatrix TEqn
(
fvm::ddt(T)
+ fvm::div(phi, T)
+ U.component(vector::X)*gamma
- fvm::laplacian(DT, T)
);

TEqn.solve();

I use cyclic boundaries for inlet and outlet and keep the mass flow constant using fvOptions.

Best regards,
Jakob

[hrvig]

Okay, nailed it..

The above equation is correct.
The values imposed at the boundaries do not matter as we are only solving for the periodic part. Afterwards the linear temperature gradient can be added and the correct field obtained. As the transport and thermodynamic properties are assumed to be independent of temperature, we are not interested in absolute temperatures.

If anyone faces the same problem for constant heat flux, please let me know and I will be more than happy to help out :-)

[dvcauwe]

Originally Posted by hrvig

Okay, nailed it..

The above equation is correct.
The values imposed at the boundaries do not matter as we are only solving for the periodic part. Afterwards the linear temperature gradient can be added and the correct field obtained. As the transport and thermodynamic properties are assumed to be independent of temperature, we are not interested in absolute temperatures.

If anyone faces the same problem for constant heat flux, please let me know and I will be more than happy to help out :-)

Exactly, once you have your periodic temperature field, the only absolute value that matters is the dimensionless temperature T+ = (Twall-T)/Ttau.

As for your previous question, yes it's normal that the bulk temperature can deviate slightly throughout the simulation due to discretization errors or so, which can be annoying when you want to gather statistics. You can prevent this by either explicitly enforcing the bulk temperature at a certain value or (I think better) by dynamically adjusting gamma just like you do your momentum source.

[David_010]

Originally Posted by hrvig

Okay, nailed it..

The above equation is correct.
The values imposed at the boundaries do not matter as we are only solving for the periodic part. Afterwards the linear temperature gradient can be added and the correct field obtained. As the transport and thermodynamic properties are assumed to be independent of temperature, we are not interested in absolute temperatures.

If anyone faces the same problem for constant heat flux, please let me know and I will be more than happy to help out :-)

Hi Jakob,

I actually need to do something similar for a pipe with periodically inserted elements and constant heat flux at the wall. Do you have any idea on how to face it?
I was thinking to add a source term to the temperature equation as you did, corresponding to the temperature gradient along the tube axis for a period. Then, iterate that term till the gradientT*cp*massFlow value is equal to the product of heat flux at the wall multiplied the tube surface -in a similar way of what is done in channelFoam for the pressure gradient and the velocity. Not sure if the idea is good.
Bests,
David

[hrvig]

Hi David,

I would start out creating a high quality mesh with inlet and outlet faces matching exactly. I ended up using blockMesh instead of snappyHexMesh even though it normally takes a little longer to create the mesh.

As dvcauwe suggested, you will end up with a bulk mean temperature that keeps increasing throughout the simulation, so I will suggest that you either explicitly force the bulk temperature to a certain value or adjust gamma throughout the simulation.

Best,
Jakob

[hcl734]

Could you upload your solver?

I am trying to modify buoyantBoussinesqSimpleFoam.

When I try to use a scalar for gamma in the solver definition like

{
    alphat = turbulence->nut()/Prt;
    alphat.correctBoundaryConditions();
    
    volScalarField alphaEff("alphaEff", turbulence->nu()/Pr + alphat);

fvScalarMatrix TEqn
(
fvm::ddt(T)
+ fvm::div(phi, T)
+ U.component(vector::X)*0.0078
      - fvm::laplacian(alphaEff, T)
     ==
        fvOptions(T)
);

    TEqn.relax();


    fvOptions.constrain(TEqn);

    TEqn.solve();

    fvOptions.correct(T);

    rhok = 1.0 - beta*(T - TRef);
}

It compiles but there is a dimension warning when starting the solver.

Using only gamma doesn't compiles, I guess because I have to define gamma somewhere but I got no clue where, cause I am pretty much a newbie to OF programming.

[hcl734]

Ok solved that on my own.

Just add a gamma entry to readtransportProperties.h

[Henry_GAO]

Hi foamers,


I'm working on a developed heat transfer (const. heat flux b.c.) with chtMultiRegionFoam.
Why above threads abandoned the jump boundary condition and turned to rewrite a solver?

I set the inlet and outlet as cyclic boundary type, and employ 'cyclic' b.c. for velocity, 'cyclicJump` for p_rgh and T.
The flow is solved first with correct p jump value.
When the flow is converged, I frozen the flow, and set the correct jump value for T, as heatFlux * Area / (phiFluid * cpFluid).
Until now, it works well.

Is there something wrong?

[Cenk20]

Hi Henry_GAO,


I'm trying almost the same simulation with buoyantPimpleFoam, but with the cyclic BC the simulation crashes. Did you have also the same problems in the beginning? How did you define the cyclic boundary conditions for p, p_rgh and T? can you maybe share your files?


Best,
Cenk

[Henry_GAO]

Hi Cenk
My simulation didn't crash.

I set the jump b.c. for p_rgh and T, And cyclic boundary for U.

The following lines are copied from 0/p_rgh



/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0/Fluid";
object p_rgh;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions [ 1 -1 -2 0 0 0 0 ];

internalField uniform 0;



boundaryField
{
...

inlet
{
type fixedJump;
patchType cyclic;
jump uniform -1;
}
outlet
{
type fixedJump;
patchType cyclic;
}
...
}

[mm66]

Originally Posted by Henry_GAO

Hi Cenk
My simulation didn't crash.

I set the jump b.c. for p_rgh and T, And cyclic boundary for U.

The following lines are copied from 0/p_rgh



/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0/Fluid";
object p_rgh;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions [ 1 -1 -2 0 0 0 0 ];

internalField uniform 0;



boundaryField
{
...

inlet
{
type fixedJump;
patchType cyclic;
jump uniform -1;
}
outlet
{
type fixedJump;
patchType cyclic;
}
...
}

Hi Henri,

May I ask how you determined the uniform jump values for T and p_rgh? I cannot understand the value of -1 in your code:

        jump            uniform -1;

Besides, I am facing some difficulty determining the jump for T in a geometry with two fluids in counter-current contact. If I am not mistaken, it should be done iteratively (similar to heat exchangers) as h is unknown in my case??