Implementation of turbulent algebraic heat flux model

Fabio1893 · January 8, 2021, 07:34

Hello dear community.

I want to implement an algebraic heat flux model, based on the one described in the StarCCM+ user guide, called Temperature Flux Model here, in OpenFOAM 8. The model is based on eddy viscosity models and needs to be implemented in a low-Re model. However, I start with an implementation in the standard k-epsilon model of OpenFOAM 8.

In OpenFOAM, if I understand correct, the heat flux is calculated, based on Boussinesq approximation, by:

$\overline{q}=- \frac{\mu_t}{Pr_t}\nabla h$ (1)

In comparison for the Temperature Flux Model, the heat flux is defined as:

$\overline{q}=-\kappa \nabla \overline{T} - \rho C_p \overline{v \theta}$ (2)

with $\overline{v\theta}$ being the turbulent heat flux, which is given as:

$\overline{v\theta}=-C_{tu0}\frac{k}{\epsilon}\left(C_{tu2}R *\nabla \overline{T}+C_{tu2} \nabla \overline{v} * \overline{v \theta}+C_{tu3}\beta\overline{\theta ^2}g \right)+C_{tu4}\left(\frac{1}{k}R-\frac{2}{3}I\right)*\overline{v\theta}$ (3)

In order to calculate the turbulent heat flux $\overline{v\theta}$ , a third transport equation for the temperature variance $\theta ^2$ needs to be solved. The transport equation for the temperature variance is very similar to the transport equation of the turbulent kinetic energy k and is defined as:

$\frac{\partial}{\partial t}\left(\rho \theta ^2\right)+\nabla * \left(\rho\theta^2\overline{v}\right)=\nabla * \left[\left(\mu+\frac{\mu_t}{\sigma_{\theta^2}}\right)\nabla\overline{\theta^2}\right]+G_\theta-\rho\epsilon_\theta$ (4)

with the production term of the temperature variance:

$G_\theta=-\rho\overline{v\theta}*\nabla\overline{T}$ (5)

Furthermore, the temperature variance dissipation rate $\epsilon_\theta$ is defined via the thermal time scale $T_\theta$

$T_\theta=\frac{\overline{\theta^2}}{2\epsilon_\theta}$ (6)

and the assumption of a constant turbulent-to-thermal time-scale ratio:

$R_T=\frac{T_\theta}{T_t}$ (7)

I have almost finished the implementation of the transport equation (4). Only the production therm (5) of the temperature variance is missing:

Code:

// Temperature variance equation
    tmp<fvScalarMatrix> theta2Eqn
    (
        fvm::ddt(alpha, rho, theta2_)
      + fvm::div(alphaRhoPhi, theta2_)
      - fvm::laplacian(alpha*rho*Dtheta2Eff(), theta2_)
      ==
          // production of theta2 is missing
      - (alpha()*rho*theta2_)/(2*R_*turbTimeScale()) // epsilonTheta2
      + theta2Source()
      + fvOptions(alpha, rho, theta2_)
    );

My question now is quite simple: How can I get access to the temperature T in my turbulence model?

EDIT:

Access to T:

Code:

const volScalarField& T_ = this->mesh_.objectRegistry::template lookupObject<volScalarField>("T");

Furthermore, I am not sure, how I could implement equation (3). I have already defined model constants, but how can I get the Reynolds stress tensor and the thermal expansion coefficient here?

If any further information is needed, please let me know. Thanks a lot in advance!

jherb · January 11, 2021, 09:51

Dear Fabio,

I implemented exactly this model in OpenFOAM (7). A German description can be found at https://www.grs.de/sites/default/files/pdf/grs-553.pdf, https://www.grs.de/publikationen/grs-553

Regarding your question about your equation (3). It is the "algebraic" equation which needs to be solved as matrix equation. So you have to re-arrange it into A x = b.
As code, this looks like:

Code:

thetaU_ = inv(leftSide) & rightSide;

Best regards,
Joachim
P.S. Are you somehow working together with NRG?

Fabio1893 · January 11, 2021, 12:02

Quote:

Originally Posted by jherb

Dear Fabio,

I implemented exactly this model in OpenFOAM (7). A German description can be found at https://www.grs.de/sites/default/files/pdf/grs-553.pdf, https://www.grs.de/publikationen/grs-553

Regarding your question about your equation (3). It is the "algebraic" equation which needs to be solved as matrix equation. So you have to re-arrange it into A x = b.
As code, this looks like:

Code:

thetaU_ = inv(leftSide) & rightSide;

Best regards,
Joachim
P.S. Are you somehow working together with NRG?

Dear Joachim,

thank you for your reply!

First of all, I want to mention that I have implemented now the thermal expansion coefficient beta, which is looked up and interpolated from a table, dependent on the temperature, which is important for my case, since I want to simulate heat transfer close to the critical point.

I try to solve the equation (3) like this:

Code:

// Turbulent heat flux equation
    tmp<fvScalarMatrix> thfEqn
    (
        thf_
      ==
      -    Ctu0_*(k()/epsilon())*(Ctu1_*this->sigma()*fvc::grad(T_) + Ctu2_*dev(twoSymm(fvc::grad(U)))*thf()+Ctu3_*beta_*theta2()*g)
      //+ Ctu4_*RSTanisotropy*thf()    
    );

    thfEqn.ref().relax();
    fvOptions.constrain(thfEqn.ref());
     solve(thfEqn);

here, sigma is the Reynolds stress tensor. However, I can not compile this code, since the types (volSymmTensorField for the therm with the Reynolds stress tensor and the term with the velocity gradient; volVectorField for the buoyancy term). The turbulent heat flux thf_ should be volScalarField, should'nt it? The anisotropy term can be neglected in my case.

So re-arranging the equation, like you have written, would solve the problem, but still I do not totally understand how to do that. Could you explain, how to re-arrange the terms into A x = b a little more detailed please?

Furthermore, I have a second question. I have implemented the transport equation for the temperature variance like that:

Code:

// Temperature variance equation
    tmp<fvScalarMatrix> theta2Eqn
    (
        fvm::ddt(alpha, rho, theta2_)
      + fvm::div(alphaRhoPhi, theta2_)
      - fvm::laplacian(alpha*rho*Dtheta2Eff(), theta2_)
      ==
        -1*alpha()*rho*thf_*mag(fvc::grad(T_)) // production of theta2
      - (alpha()*rho*theta2_)/(2*R_*turbTimeScale()) // epsilonTheta2
      + theta2Source()
      + fvOptions(alpha, rho, theta2_)
    );

Here, first, I tried to write the production term of the temperature variance as

Code:

-1*alpha()*rho*thf_*fvc::grad(T_)

instead of

Code:

-1*alpha()*rho*thf_*mag(fvc::grad(T_))

but I got an error due to the types. So the question is now, is it correct to use the mag(), to yield a volScalarField type for the production term?

EDIT (SOLVED):

of course, thf_ needs to be a volVectorField. Therefore the temperature variance can be implemented as follows:

Code:

// Temperature variance equation
    tmp<fvScalarMatrix> theta2Eqn
    (
        fvm::ddt(alpha, rho, theta2_)
      + fvm::div(alphaRhoPhi, theta2_)
      + alpha()*rho*thf_&fvc::gradT_) // production of theta2
      + alpha()*(rho*theta2_)/(2*R_*turbTimeScale()) // epsilonTheta2
      ==
        fvm::laplacian(alpha*rho*Dtheta2Eff(), theta2_)
      + theta2Source()
      + fvOptions(alpha, rho, theta2_)
    );

PS: No, I have nothing to do with GRS, I am PhD student, working on SCWR.

Thanks a lot so far!

Fabio1893 · January 12, 2021, 09:59

I am still struggling with the algebraic equation. I wrote it as vector equation, but get a type error while compiling it, when I do it like this:

Code:

//Turbulen heat flux equation
    tmp<fvVectorMatrix> thfEqn
    (
        thf_
      ==
      - Ctu0_*(k()/epsilon())*(Ctu1_*this->sigma() & fvc::grad(T_)
      + Ctu2_*dev(twoSymm(fvc::grad(U))) & thf()
      + Ctu3_*beta_*theta2()*g)
      + Ctu4_*RSTanisotropy & thf()
    );

However, it is possible to define all the terms as volVectorFields like this:

Code:

volScalarField C = Ctu0_*(k()/epsilon());

    volVectorField Term1 = C*Ctu1_*this->sigma() & fvc::grad(T_);
    volVectorField Term2 = C*Ctu2_*dev(twoSymm(fvc::grad(U))) & thf();
    volVectorField Term3 = C*Ctu3_*beta_*theta2()*g;
    volVectorField Term4 = Ctu4_*RSTanisotropy & thf();

Now, it should be possible to write the equation without type errors:

Code:

//Turbulent heat flux equation
    tmp<fvVectorMatrix> thfEqn
    (
        thf_
      ==
      - Term1
      - Term2
      - Term3 
      + Term4
    );

    thfEqn.ref().relax();
    fvOptions.constrain(thfEqn.ref());
    solve(thfEqn);

However, I get the following error:

Code:

error: invalid use of void expression
  481 |     tmp<fvVectorMatrix> thfEqn

Does someone have any idea?

jherb · January 16, 2021, 14:08

First to the problem of the algebraic equation.

$\overline{v\theta}=-C_{tu0}\frac{k}{\epsilon}\left(C_{tu2}R *\nabla \overline{T}+C_{tu2} \nabla \overline{v} * \overline{v \theta}+C_{tu3}\beta\overline{\theta ^2}g \right)+C_{tu4}\left(\frac{1}{k}R-\frac{2}{3}I\right)*\overline{v\theta}$ (3)

This is the algebraic equation. You have to rearrange it, so that the matrix equation is clear:

$\left(I + C_{tu0}\frac{k}{\epsilon} C_{tu2} \nabla \overline{v} - C_{tu4}\left(\frac{1}{k}R-\frac{2}{3}I\right) \right) * \overline{v\theta}=-C_{tu0}\frac{k}{\epsilon}\left(C_{tu2}R *\nabla \overline{T}+C_{tu3}\beta\overline{\theta ^2}g \right)$

So the matrix A is everything before $\overline{v\theta}$ . It is of type volTensorField. So you can assign it to such a variable, e. g.:

Code:

    tmp<volTensorField> tgradU2 = fvc::grad(U_);
    const volTensorField& gradU2 = tgradU2();


    volTensorField leftSide
    (
        I
      + Ctu0_*k_/epsilon_*Ctu2_*gradU2
      - Ct4_*...
      );

The right hand side is a volVectorField:

Code:

    volVectorField gradT(fvc::grad(T));

    volVectorField rightSide
    (
        -Ctu0_*k_/epsilon_*
        (
            (Ctu2_*R & gradT)
          +
            (Ctu3_*beta_*g*thetaSqr_)
        )
    );

The matrix equation can then be solved for $\overline{v\theta}$ :

Code:

thetaU_ = inv(leftSide) & rightSide;

By the way, thetaU_ is also a volVectorField.

jherb · January 16, 2021, 14:25

Which solver will you use? Are you aware, that you have to modify the energy/temperature equation of the solver? Now the heat/energy flux becomes becomes a sum of the convective heat flux and the turbulent one, calculated by thetaU_. See e. g. equation 2 in https://www.researchgate.net/publica...PRANDTL_FLUIDS

Fabio1893 · January 19, 2021, 05:37

Thank you a lot again, the model compiles and runs now. - A first try, for now, without any use, since I do not use the thetaU yet. -

I use a slightly modified solver based on rhoSimpleFoam. I have added several fields in which I am interested and for simulation of mixed convection cases, I found that an additional buoyancy term in the momentum equation increases the simulation accuracy. For forced convection cases however, I use the "standard" momentum equation of rhoSimpleFoam.

My idea was to redefine the calculation of the heat flux in the eddyDiffusivity.C, using thetaU, instead of changing the solvers energy equation. I think this will not make any difference at all, but maybe you have experience with that? I try to go both ways and will report it here, if someone is interested.

jherb · January 19, 2021, 16:11

I also added buoyancy terms to the turbulence equations. And for the other fields it is included in buoyantBoussineqPimpleFoam.

I made the experience that the simulations are unstable at the beginning and I had to use some tricks like limiting terms. When the simulations reached a steady state, these tricks were not longer needed. I think the problem is some terms on the right hand side of equation (4) are transport terms but do not have the "standard" form $\left(\rho\overline{v}\psi\right)$ or $\left(\phi\psi\right)$ .

Fabio1893 · January 20, 2021, 04:42

I have also added buoyancy terms for k and epsilon, like this:

$G_{bk} = -\rho\beta g \overline{v \theta}$

and

$G_{b \epsilon} = -\left(\frac{\epsilon}{k}\right)C_{\epsilon}\rho\beta g \overline{v \theta}$

I have not implemented the last "step" yet, to use thetaU in the energy equation or in the eddy diffusivity model. Should be quite easy, but I did not find out yet, how to access the thetaU field from outside the library of the turbulence model. Anyway thanks for the hint with limiting terms.

What boundary conditions do you use for the thetaSqr?

Andrea1984 · January 22, 2021, 05:13

Hi Joachim,

do you have any plans to release the source code of your implementation? Google translate is making a terrible job with your document

Regards,
Andrea

Quote:

Originally Posted by jherb

Dear Fabio,

I implemented exactly this model in OpenFOAM (7). A German description can be found at https://www.grs.de/sites/default/files/pdf/grs-553.pdf, https://www.grs.de/publikationen/grs-553

Regarding your question about your equation (3). It is the "algebraic" equation which needs to be solved as matrix equation. So you have to re-arrange it into A x = b.
As code, this looks like:

Code:

thetaU_ = inv(leftSide) & rightSide;

Best regards,
Joachim
P.S. Are you somehow working together with NRG?

jherb · January 22, 2021, 16:38

Currently, the code is not really good enough to open source it. And it might get very complicated to get this aproved by the management.

But I might be able to answer specific implementation questions here.

Andrea1984 · January 23, 2021, 12:20

No problem, I fully understand your point.

BTW I have been working at NRG on this stuff for a while

Fabio1893 · January 26, 2021, 07:23

Quote:

Originally Posted by jherb

Which solver will you use? Are you aware, that you have to modify the energy/temperature equation of the solver? Now the heat/energy flux becomes becomes a sum of the convective heat flux and the turbulent one, calculated by thetaU_. See e. g. equation 2 in https://www.researchgate.net/publica...PRANDTL_FLUIDS

Hi Joachim,

since the energy equation is defined as fvScalarMatrix and thetaU_ is a volVectorField, I still can not finish my implementation. The same problem appears if I try to modify alphat in the eddyDiffusivity.C, which is a volScalarField. Can you give me any hint how to solve this problem?

Best regards,
Fabian

jherb · January 29, 2021, 18:50

This is the energy equation you want to change?
https://github.com/OpenFOAM/OpenFOAM...oam/EEqn.H#L12

In the new form introduced in OpenFOAM-8 the diffusion is implemented in this divq method. You have to search for it. I am not sure, but is probably this implementation you need:
https://github.com/OpenFOAM/OpenFOAM.../Fourier.C#L88
(only one fluid component)

Here you see the use of the Reynolds analogy. Now you have to replace this term in the energy equation with the new diffusion term, e.g. the second summand on the right hand side of eq. (2) in https://www.researchgate.net/publica...PRANDTL_FLUIDS

Quote:

Originally Posted by Fabio1893

Hi Joachim,

since the energy equation is defined as fvScalarMatrix and thetaU_ is a volVectorField, I still can not finish my implementation. The same problem appears if I try to modify alphat in the eddyDiffusivity.C, which is a volScalarField. Can you give me any hint how to solve this problem?

Best regards,
Fabian

Fabio1893 · February 2, 2021, 09:53

Hi Joachim,

yes, this energy equation is exactly the one I want to modify, with my changes marked in red.

Code:

{
     volScalarField& he = thermo.he();
     volVectorField thetaU = U.db().objectRegistry::template lookupObject<volVectorField>("thetaU");
 
    fvScalarMatrix EEqn
    (
        fvm::div(phi, he)
      + (
            he.name() == "e"
          ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
          : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
        )
       //+ thermophysicalTransport->divq(he)
       - rho*thermo.Cp()*fvc::div(thetaU)
      ==
        fvOptions(rho, he)
    );

    EEqn.relax();

    fvOptions.constrain(EEqn);

    EEqn.solve();

    fvOptions.correct(he);

    thermo.correct();
}

The divq(he) function is defined in eddyDiffusivity.C:

Code:

template<class TurbulenceThermophysicalTransportModel>
tmp<fvScalarMatrix>
eddyDiffusivity<TurbulenceThermophysicalTransportModel>::divq
(
    volScalarField& he
) const
{
    return -fvm::laplacian(this->alpha()*this->alphaEff(), he);
 }

However, my implenetation seems to be incorrect, since I muss the first term of equation (2). Therefore I think, alpha()*alphaEff() represents the part in brackets of this equation:
$q=-\left(\kappa+\frac{\mu_t C_p}{Pr_t}\right)\nabla T$
since in the algebraic model, the left part is still multiplied by the temperature gradient and the right part is something multiplied by thetaU (see eq (2)), I have to "split" this alpha()*alphaEff() part.

Best regards,
Fabian

jherb · February 2, 2021, 19:07

I acctually used buoyantBoussinesqPimpleFoam which has a temperature equation. So I could directly use equation (2).

https://github.com/OpenFOAM/OpenFOAM...pleFoam/TEqn.H

There

Code:

fvm::laplacian(alphaEff, T)

becomes something like

Code:

fvm::laplacian(turbulence->nu()/Pr, T) + fvc::div(thetaU)

Now you have to somehow re-write this for the energy equation. A problem might be, that alpha/cp/... might not be constant over space, so you have to apply the chain rule. But I have not thought about this yet.

jherb · February 2, 2021, 19:28

Looking at equation (3) of https://www.sciencedirect.com/scienc...38573320309931
and considering, that https://github.com/OpenFOAM/OpenFOAM.../Fourier.C#L88 is just the laminar part the diffusion term should now be something like:

Code:

fvm::laplacian(this->alpha()*this->thermo().alpha(), he) + rho*thermo.Cp()*fvc::div(thetaU)

(I do not understand, from were the factor this->alpha()*this->thermo().alpha() comes)

Edit: I think the first alpha is the phase fraction (because the model can also used for multi-phase simulations) and the second alpha (thermo().alpha()) is the thermal diffusivity. So for single phase flows, the first alpha = 1

Fabio1893 · February 4, 2021, 07:06

Hi Joachim,

yes you are right, the first alpha is the phase fraction, this alpha appears also in all of the transport equations of the turbulence model. The second alpha, which is called alphaEff in the divq() function in my previous post, is, as you said, thermal diffusivity. We can find the definition in the heThermo.C file:

Code:

template<class BasicThermo, class MixtureType>
Foam::tmp<Foam::volScalarField>
Foam::heThermo<BasicThermo, MixtureType>::alphaEff
(
    const volScalarField& alphat
) const
{
    return volScalarField::New
    (
        "alphaEff",
        this->CpByCpv()*(this->alpha_ + alphat)
    );
}

in here, the turbulent part alphat comes from eddyDiffusivity.C:

Code:

template<class TurbulenceThermophysicalTransportModel>
void eddyDiffusivity<TurbulenceThermophysicalTransportModel>::correctAlphat()
{
    alphat_ =
        this->momentumTransport().rho()
       *this->momentumTransport().nut()/Prt_;
    alphat_.correctBoundaryConditions();
}

and alpha is the laminar part. CpByCpv equals 1 on my case (I am quite sure at least), since I use sensible enthalpy model.

Therefore, similar to your suggestion, I write the EEqn like this:

Code:

fvScalarMatrix EEqn
    (
        fvm::div(phi, he)
      + (
            he.name() == "e"
          ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
          : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
        )
     + thermophysicalTransport->divq(he) - rho*thermo.Cp()*fvc::div(thetaU)
     ==
        fvOptions(rho, he)
    );

So the laplacian is in the divq function and I add the "algebraic" part.

However, rho and Cp are definitely not constant over space, since I have strong variations in this properties with temperature close to the critical point and I have quite high spatial temperature gradients in my simulation cases.

For now, the code is not running stable. I always get bounding thetaSqr, continuity errors explode after a while. Relaxation factors and numerical schemes have some influence, but I did not have a real success yet. I have tested the same case with the same turbulence model, without the additional transport equation and with the "standard" EEqn in the solver. Runs without any problem and results are what one would expect.

Another point is, that the code only runs on a single core, but not in parallel. This seems to be due to the change in the energy equation.

Best regards,
Fabian

jherb · February 7, 2021, 12:04

I think, this is not correct:

Code:

thermophysicalTransport->divq(he) - rho*thermo.Cp()*fvc::div(thetaU)

because the divq part also includes the alphat part. But his part should be replace by the thetaU...

Regarding the instabilities: This is exactly, what I saw. I think the problem might be related to the additional source terms in the different turbulence equations. Those terms are physically some kind of transport terms, but they are not in the form of fvc::div(phi, XXX). If they were one could use upstream schemes. But so, it looks like more of quantity of a cell is removed, compared to what is stored in the cell.

Fabio1893 · February 11, 2021, 10:24

Yes you are right about that, since alphaEff is the molecular and the turbulent part together.

Therefore, I define a new function "divAlgebraic", which replaces the divq function in the EEqn of the solver. Instead of alphaEff, alphahe is used here, which is only the molecular part. Furthermore, I add the rho*Cp*div(thetaU) for the algebraic heat flux model, and also add the alpha for the phase fraction.

Code:

template<class TurbulenceThermophysicalTransportModel>
tmp<fvScalarMatrix>
eddyDiffusivity<TurbulenceThermophysicalTransportModel>::divAlgebraic
(
    volScalarField& he,
    volVectorField& thetaU
) const
{
    return -fvm::laplacian(this->alpha()*this->thermo().alphahe(), he) - this->alpha()*this->momentumTransport().rho()*this->thermo().Cp()*fvc::div(thetaU);
}

The definition of alphahe comes from heThermo.C:

Code:

template<class BasicThermo, class MixtureType>
Foam::tmp<Foam::volScalarField>
Foam::heThermo<BasicThermo, MixtureType>::alphahe() const
{
    return volScalarField::New
    (
        "alphahe",
        this->CpByCpv()*this->alpha_
    );
 }

I think now the equations are correct, but still it is only possible to run on one single core, so there must be anything going wrong in the implementation.

The biggest problem however, is the stability of the simulation. I understand your point about the form of the transport terms, but I dont know how to solve the issue. Can you tell me which divergence schemes you use?

Best regads,
Fabian

January 8, 2021, 07:34	Implementation of turbulent algebraic heat flux model	#1
Fabio1893 New Member Join Date: Jan 2021 Posts: 10 Rep Power: 5	Hello dear community. I want to implement an algebraic heat flux model, based on the one described in the StarCCM+ user guide, called Temperature Flux Model here, in OpenFOAM 8. The model is based on eddy viscosity models and needs to be implemented in a low-Re model. However, I start with an implementation in the standard k-epsilon model of OpenFOAM 8. In OpenFOAM, if I understand correct, the heat flux is calculated, based on Boussinesq approximation, by: $\overline{q}=- \frac{\mu_t}{Pr_t}\nabla h$ (1) In comparison for the Temperature Flux Model, the heat flux is defined as: $\overline{q}=-\kappa \nabla \overline{T} - \rho C_p \overline{v \theta}$ (2) with $\overline{v\theta}$ being the turbulent heat flux, which is given as: $\overline{v\theta}=-C_{tu0}\frac{k}{\epsilon}\left(C_{tu2}R \nabla \overline{T}+C_{tu2} \nabla \overline{v} \overline{v \theta}+C_{tu3}\beta\overline{\theta ^2}g \right)+C_{tu4}\left(\frac{1}{k}R-\frac{2}{3}I\right)\overline{v\theta}$ (3) In order to calculate the turbulent heat flux $\overline{v\theta}$ , a third transport equation for the temperature variance $\theta ^2$ needs to be solved. The transport equation for the temperature variance is very similar to the transport equation of the turbulent kinetic energy k and is defined as: $\frac{\partial}{\partial t}\left(\rho \theta ^2\right)+\nabla \left(\rho\theta^2\overline{v}\right)=\nabla * \left[\left(\mu+\frac{\mu_t}{\sigma_{\theta^2}}\right)\nabla\overline{\theta^2}\right]+G_\theta-\rho\epsilon_\theta$ (4) with the production term of the temperature variance: $G_\theta=-\rho\overline{v\theta}\nabla\overline{T}$ (5) Furthermore, the temperature variance dissipation rate $\epsilon_\theta$ is defined via the thermal time scale $T_\theta$ $T_\theta=\frac{\overline{\theta^2}}{2\epsilon_\theta}$ (6) and the assumption of a constant turbulent-to-thermal time-scale ratio: $R_T=\frac{T_\theta}{T_t}$ (7) I have almost finished the implementation of the transport equation (4). Only the production therm (5) of the temperature variance is missing: Code: // Temperature variance equation tmp<fvScalarMatrix> theta2Eqn ( fvm::ddt(alpha, rho, theta2_) + fvm::div(alphaRhoPhi, theta2_) - fvm::laplacian(alpharhoDtheta2Eff(), theta2_) == // production of theta2 is missing - (alpha()rhotheta2_)/(2R_turbTimeScale()) // epsilonTheta2 + theta2Source() + fvOptions(alpha, rho, theta2_) ); My question now is quite simple: How can I get access to the temperature T in my turbulence model? EDIT: Access to T: Code: const volScalarField& T_ = this->mesh_.objectRegistry::template lookupObject<volScalarField>("T"); Furthermore, I am not sure, how I could implement equation (3). I have already defined model constants, but how can I get the Reynolds stress tensor and the thermal expansion coefficient here? If any further information is needed, please let me know. Thanks a lot in advance! Minghao_Li and lpz456 like this. Last edited by Fabio1893; January 8, 2021 at 08:38.*

January 11, 2021, 09:51		#2
jherb Senior Member Joachim Herb Join Date: Sep 2010 Posts: 650 Rep Power: 22	Dear Fabio, I implemented exactly this model in OpenFOAM (7). A German description can be found at https://www.grs.de/sites/default/files/pdf/grs-553.pdf, https://www.grs.de/publikationen/grs-553 Regarding your question about your equation (3). It is the "algebraic" equation which needs to be solved as matrix equation. So you have to re-arrange it into A x = b. As code, this looks like: Code: thetaU_ = inv(leftSide) & rightSide; Best regards, Joachim P.S. Are you somehow working together with NRG?

January 16, 2021, 14:08		#5
jherb Senior Member Joachim Herb Join Date: Sep 2010 Posts: 650 Rep Power: 22	First to the problem of the algebraic equation. $\overline{v\theta}=-C_{tu0}\frac{k}{\epsilon}\left(C_{tu2}R \nabla \overline{T}+C_{tu2} \nabla \overline{v} \overline{v \theta}+C_{tu3}\beta\overline{\theta ^2}g \right)+C_{tu4}\left(\frac{1}{k}R-\frac{2}{3}I\right)\overline{v\theta}$ (3) This is the algebraic equation. You have to rearrange it, so that the matrix equation is clear: $\left(I + C_{tu0}\frac{k}{\epsilon} C_{tu2} \nabla \overline{v} - C_{tu4}\left(\frac{1}{k}R-\frac{2}{3}I\right) \right) \overline{v\theta}=-C_{tu0}\frac{k}{\epsilon}\left(C_{tu2}R \nabla \overline{T}+C_{tu3}\beta\overline{\theta ^2}g \right)$ So the matrix A is everything before $\overline{v\theta}$ . It is of type volTensorField. So you can assign it to such a variable, e. g.: Code: tmp<volTensorField> tgradU2 = fvc::grad(U_); const volTensorField& gradU2 = tgradU2(); volTensorField leftSide ( I + Ctu0_k_/epsilon_Ctu2_gradU2 - Ct4_... ); The right hand side is a volVectorField: Code: volVectorField gradT(fvc::grad(T)); volVectorField rightSide ( -Ctu0_k_/epsilon_* ( (Ctu2_R & gradT) + (Ctu3_beta_gthetaSqr_) ) ); The matrix equation can then be solved for $\overline{v\theta}$ : Code: thetaU_ = inv(leftSide) & rightSide; By the way, thetaU_ is also a volVectorField. Farid and lpz456 like this.

January 19, 2021, 05:37		#7
Fabio1893 New Member Join Date: Jan 2021 Posts: 10 Rep Power: 5	Thank you a lot again, the model compiles and runs now. - A first try, for now, without any use, since I do not use the thetaU yet. - I use a slightly modified solver based on rhoSimpleFoam. I have added several fields in which I am interested and for simulation of mixed convection cases, I found that an additional buoyancy term in the momentum equation increases the simulation accuracy. For forced convection cases however, I use the "standard" momentum equation of rhoSimpleFoam. My idea was to redefine the calculation of the heat flux in the eddyDiffusivity.C, using thetaU, instead of changing the solvers energy equation. I think this will not make any difference at all, but maybe you have experience with that? I try to go both ways and will report it here, if someone is interested. lpz456 likes this.

January 19, 2021, 16:11		#8
jherb Senior Member Joachim Herb Join Date: Sep 2010 Posts: 650 Rep Power: 22	I also added buoyancy terms to the turbulence equations. And for the other fields it is included in buoyantBoussineqPimpleFoam. I made the experience that the simulations are unstable at the beginning and I had to use some tricks like limiting terms. When the simulations reached a steady state, these tricks were not longer needed. I think the problem is some terms on the right hand side of equation (4) are transport terms but do not have the "standard" form $\left(\rho\overline{v}\psi\right)$ or $\left(\phi\psi\right)$ . lpz456 likes this.

January 16, 2021, 14:25		#6
jherb Senior Member Joachim Herb Join Date: Sep 2010 Posts: 650 Rep Power: 22	Which solver will you use? Are you aware, that you have to modify the energy/temperature equation of the solver? Now the heat/energy flux becomes becomes a sum of the convective heat flux and the turbulent one, calculated by thetaU_. See e. g. equation 2 in https://www.researchgate.net/publica...PRANDTL_FLUIDS

January 20, 2021, 04:42		#9
Fabio1893 New Member Join Date: Jan 2021 Posts: 10 Rep Power: 5	I have also added buoyancy terms for k and epsilon, like this: $G_{bk} = -\rho\beta g \overline{v \theta}$ and $G_{b \epsilon} = -\left(\frac{\epsilon}{k}\right)C_{\epsilon}\rho\beta g \overline{v \theta}$ I have not implemented the last "step" yet, to use thetaU in the energy equation or in the eddy diffusivity model. Should be quite easy, but I did not find out yet, how to access the thetaU field from outside the library of the turbulence model. Anyway thanks for the hint with limiting terms. What boundary conditions do you use for the thetaSqr?

January 22, 2021, 16:38		#11
jherb Senior Member Joachim Herb Join Date: Sep 2010 Posts: 650 Rep Power: 22	Currently, the code is not really good enough to open source it. And it might get very complicated to get this aproved by the management. But I might be able to answer specific implementation questions here.

January 23, 2021, 12:20		#12
Andrea1984 Senior Member Andrea Join Date: Feb 2012 Location: Leeds, UK Posts: 179 Rep Power: 16	No problem, I fully understand your point. BTW I have been working at NRG on this stuff for a while

February 2, 2021, 09:53		#15
Fabio1893 New Member Join Date: Jan 2021 Posts: 10 Rep Power: 5	Hi Joachim, yes, this energy equation is exactly the one I want to modify, with my changes marked in red. Code: { volScalarField& he = thermo.he(); volVectorField thetaU = U.db().objectRegistry::template lookupObject<volVectorField>("thetaU"); fvScalarMatrix EEqn ( fvm::div(phi, he) + ( he.name() == "e" ? fvc::div(phi, volScalarField("Ekp", 0.5magSqr(U) + p/rho)) : fvc::div(phi, volScalarField("K", 0.5magSqr(U))) ) //+ thermophysicalTransport->divq(he) - rhothermo.Cp()fvc::div(thetaU) == fvOptions(rho, he) ); EEqn.relax(); fvOptions.constrain(EEqn); EEqn.solve(); fvOptions.correct(he); thermo.correct(); } The divq(he) function is defined in eddyDiffusivity.C: Code: template<class TurbulenceThermophysicalTransportModel> tmp<fvScalarMatrix> eddyDiffusivity<TurbulenceThermophysicalTransportModel>::divq ( volScalarField& he ) const { return -fvm::laplacian(this->alpha()this->alphaEff(), he); } However, my implenetation seems to be incorrect, since I muss the first term of equation (2). Therefore I think, alpha()alphaEff() represents the part in brackets of this equation: $q=-\left(\kappa+\frac{\mu_t C_p}{Pr_t}\right)\nabla T$ since in the algebraic model, the left part is still multiplied by the temperature gradient and the right part is something multiplied by thetaU (see eq (2)), I have to "split" this alpha()*alphaEff() part. Best regards, Fabian

February 2, 2021, 19:07		#16
jherb Senior Member Joachim Herb Join Date: Sep 2010 Posts: 650 Rep Power: 22	I acctually used buoyantBoussinesqPimpleFoam which has a temperature equation. So I could directly use equation (2). https://github.com/OpenFOAM/OpenFOAM...pleFoam/TEqn.H There Code: fvm::laplacian(alphaEff, T) becomes something like Code: fvm::laplacian(turbulence->nu()/Pr, T) + fvc::div(thetaU) Now you have to somehow re-write this for the energy equation. A problem might be, that alpha/cp/... might not be constant over space, so you have to apply the chain rule. But I have not thought about this yet.

February 2, 2021, 19:28		#17
jherb Senior Member Joachim Herb Join Date: Sep 2010 Posts: 650 Rep Power: 22	Looking at equation (3) of https://www.sciencedirect.com/scienc...38573320309931 and considering, that https://github.com/OpenFOAM/OpenFOAM.../Fourier.C#L88 is just the laminar part the diffusion term should now be something like: Code: fvm::laplacian(this->alpha()this->thermo().alpha(), he) + rhothermo.Cp()fvc::div(thetaU) (I do not understand, from were the factor this->alpha()this->thermo().alpha() comes) Edit: I think the first alpha is the phase fraction (because the model can also used for multi-phase simulations) and the second alpha (thermo().alpha()) is the thermal diffusivity. So for single phase flows, the first alpha = 1 Last edited by jherb; February 3, 2021 at 15:12.

February 4, 2021, 07:06		#18
Fabio1893 New Member Join Date: Jan 2021 Posts: 10 Rep Power: 5	Hi Joachim, yes you are right, the first alpha is the phase fraction, this alpha appears also in all of the transport equations of the turbulence model. The second alpha, which is called alphaEff in the divq() function in my previous post, is, as you said, thermal diffusivity. We can find the definition in the heThermo.C file: Code: template<class BasicThermo, class MixtureType> Foam::tmp<Foam::volScalarField> Foam::heThermo<BasicThermo, MixtureType>::alphaEff ( const volScalarField& alphat ) const { return volScalarField::New ( "alphaEff", this->CpByCpv()(this->alpha_ + alphat) ); } in here, the turbulent part alphat comes from eddyDiffusivity.C: Code: template<class TurbulenceThermophysicalTransportModel> void eddyDiffusivity<TurbulenceThermophysicalTransportModel>::correctAlphat() { alphat_ = this->momentumTransport().rho() this->momentumTransport().nut()/Prt_; alphat_.correctBoundaryConditions(); } and alpha is the laminar part. CpByCpv equals 1 on my case (I am quite sure at least), since I use sensible enthalpy model. Therefore, similar to your suggestion, I write the EEqn like this: Code: fvScalarMatrix EEqn ( fvm::div(phi, he) + ( he.name() == "e" ? fvc::div(phi, volScalarField("Ekp", 0.5magSqr(U) + p/rho)) : fvc::div(phi, volScalarField("K", 0.5magSqr(U))) ) + thermophysicalTransport->divq(he) - rhothermo.Cp()fvc::div(thetaU) == fvOptions(rho, he) ); So the laplacian is in the divq function and I add the "algebraic" part. However, rho and Cp are definitely not constant over space, since I have strong variations in this properties with temperature close to the critical point and I have quite high spatial temperature gradients in my simulation cases. For now, the code is not running stable. I always get bounding thetaSqr, continuity errors explode after a while. Relaxation factors and numerical schemes have some influence, but I did not have a real success yet. I have tested the same case with the same turbulence model, without the additional transport equation and with the "standard" EEqn in the solver. Runs without any problem and results are what one would expect. Another point is, that the code only runs on a single core, but not in parallel. This seems to be due to the change in the energy equation. Best regards, Fabian

February 7, 2021, 12:04		#19
jherb Senior Member Joachim Herb Join Date: Sep 2010 Posts: 650 Rep Power: 22	I think, this is not correct: Code: thermophysicalTransport->divq(he) - rhothermo.Cp()fvc::div(thetaU) because the divq part also includes the alphat part. But his part should be replace by the thetaU... Regarding the instabilities: This is exactly, what I saw. I think the problem might be related to the additional source terms in the different turbulence equations. Those terms are physically some kind of transport terms, but they are not in the form of fvc::div(phi, XXX). If they were one could use upstream schemes. But so, it looks like more of quantity of a cell is removed, compared to what is stored in the cell.

February 11, 2021, 10:24		#20
Fabio1893 New Member Join Date: Jan 2021 Posts: 10 Rep Power: 5	Yes you are right about that, since alphaEff is the molecular and the turbulent part together. Therefore, I define a new function "divAlgebraic", which replaces the divq function in the EEqn of the solver. Instead of alphaEff, alphahe is used here, which is only the molecular part. Furthermore, I add the rhoCpdiv(thetaU) for the algebraic heat flux model, and also add the alpha for the phase fraction. Code: template<class TurbulenceThermophysicalTransportModel> tmp<fvScalarMatrix> eddyDiffusivity<TurbulenceThermophysicalTransportModel>::divAlgebraic ( volScalarField& he, volVectorField& thetaU ) const { return -fvm::laplacian(this->alpha()this->thermo().alphahe(), he) - this->alpha()this->momentumTransport().rho()this->thermo().Cp()fvc::div(thetaU); } The definition of alphahe comes from heThermo.C: Code: template<class BasicThermo, class MixtureType> Foam::tmp<Foam::volScalarField> Foam::heThermo<BasicThermo, MixtureType>::alphahe() const { return volScalarField::New ( "alphahe", this->CpByCpv()*this->alpha_ ); } I think now the equations are correct, but still it is only possible to run on one single core, so there must be anything going wrong in the implementation. The biggest problem however, is the stability of the simulation. I understand your point about the form of the transport terms, but I dont know how to solve the issue. Can you tell me which divergence schemes you use? Best regads, Fabian lpz456 likes this.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Centrifugal fan-reverse flow in outlet lesds to a mass in flow field	xiexing	CFX	3	March 29, 2017 11:00
Wrong flow in ratating domain problem	Sanyo	CFX	17	August 15, 2015 07:20
CFX radiatin model wall heat flux imblance	cscfx	CFX	0	May 21, 2014 04:39
UDF for Heat Exchanger model	francois louw	FLUENT	2	July 16, 2010 03:21
Natural convection - Inlet boundary condition	max91	CFX	1	July 29, 2008 21:28