# Implementation of turbulent algebraic heat flux model

 Register Blogs Members List Search Today's Posts Mark Forums Read

 January 8, 2021, 06:34 Implementation of turbulent algebraic heat flux model #1 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: (1) In comparison for the Temperature Flux Model, the heat flux is defined as: (2) with being the turbulent heat flux, which is given as: (3) In order to calculate the turbulent heat flux , a third transport equation for the temperature variance 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: (4) with the production term of the temperature variance: (5) Furthermore, the temperature variance dissipation rate is defined via the thermal time scale (6) and the assumption of a constant turbulent-to-thermal time-scale ratio: (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 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("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 likes this. Last edited by Fabio1893; January 8, 2021 at 07:38.

 January 11, 2021, 08:51 #2 Senior Member   Joachim Herb Join Date: Sep 2010 Posts: 650 Rep Power: 21 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 11, 2021, 11:02
#3
New Member

Join Date: Jan 2021
Posts: 10
Rep Power: 5
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,

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_
==
//+ 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_)
==
- (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_)
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!

Last edited by Fabio1893; January 12, 2021 at 08:54.

 January 12, 2021, 08:59 #4 New Member   Join Date: Jan 2021 Posts: 10 Rep Power: 5 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 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 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 thfEqn Does someone have any idea? Last edited by Fabio1893; January 12, 2021 at 10:05.

 January 16, 2021, 13:08 #5 Senior Member   Joachim Herb Join Date: Sep 2010 Posts: 650 Rep Power: 21 First to the problem of the algebraic equation. (3) This is the algebraic equation. You have to rearrange it, so that the matrix equation is clear: So the matrix A is everything before . It is of type volTensorField. So you can assign it to such a variable, e. g.: Code:  tmp 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 : Code: thetaU_ = inv(leftSide) & rightSide; By the way, thetaU_ is also a volVectorField. Farid likes this.

 January 16, 2021, 13:25 #6 Senior Member   Joachim Herb Join Date: Sep 2010 Posts: 650 Rep Power: 21 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 19, 2021, 04:37 #7 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.

 January 19, 2021, 15:11 #8 Senior Member   Joachim Herb Join Date: Sep 2010 Posts: 650 Rep Power: 21 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 or .

 January 20, 2021, 03:42 #9 New Member   Join Date: Jan 2021 Posts: 10 Rep Power: 5 I have also added buoyancy terms for k and epsilon, like this: and 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, 04:13
#10
Senior Member

Andrea
Join Date: Feb 2012
Location: Leeds, UK
Posts: 179
Rep Power: 16
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?

 January 22, 2021, 15:38 #11 Senior Member   Joachim Herb Join Date: Sep 2010 Posts: 650 Rep Power: 21 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, 11:20 #12 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

January 26, 2021, 06:23
#13
New Member

Join Date: Jan 2021
Posts: 10
Rep Power: 5
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

January 29, 2021, 17:50
#14
Senior Member

Joachim Herb
Join Date: Sep 2010
Posts: 650
Rep Power: 21
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

 February 2, 2021, 08:53 #15 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("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 tmp eddyDiffusivity::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: 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, 18:07 #16 Senior Member   Joachim Herb Join Date: Sep 2010 Posts: 650 Rep Power: 21 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, 18:28 #17 Senior Member   Joachim Herb Join Date: Sep 2010 Posts: 650 Rep Power: 21 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 Last edited by jherb; February 3, 2021 at 14:12.

 February 4, 2021, 06:06 #18 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 Foam::tmp Foam::heThermo::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 void eddyDiffusivity::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

 February 7, 2021, 11:04 #19 Senior Member   Joachim Herb Join Date: Sep 2010 Posts: 650 Rep Power: 21 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.

 February 11, 2021, 09:24 #20 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 rho*Cp*div(thetaU) for the algebraic heat flux model, and also add the alpha for the phase fraction. Code: template tmp eddyDiffusivity::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 Foam::tmp Foam::heThermo::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