Implementation of turbulent algebraic heat flux model

agustinvo · February 15, 2021, 14:03

Hi @Andrea, nice to read you around!

I managed to implement this some months ago in OpenFOAM 7, and it provides reasonable results in forced convection. As they say before, it seems quite unstable, so it is better to initialize with another turbulence model. If you are trying to use it in natural convection, it can be painful.

About the fvSolutions, I use to relax the THF (check any pEqn.H to check out how you can do this), and try to use upwind schemes, at least at the beginning of the simulation.

The problems I found out with this model is the numericla stability: the leftSide as @jherb refers to can have a very bad condition number, which leads to spurious THF values, even if the region is isothermal. I'm looking for some workarounds for this.

Good luck!

BTW: what is the test case you use to validate your implementation?

Andrea1984 · February 15, 2021, 14:15

If it can be of any consolation, this turbulent heat flux closure was pretty unstable in both STAR-CCM+ and Code_Saturne as well.

On top of what Augustin mentioned above, I would also suggest to ramp-up the values of the coefficients gradually from ~ 0 to their final value. It might also be helpful to set Ct4 = 0 altogether (unless the term multiplied by Ct4 plays a key role in the physics of the case that you are modelling).

Good luck,
Andrea

jherb · February 16, 2021, 16:05

This exactly what I also saw: At the beginning of the simulation of a Rayleigh–Bénard convection, when the whole flow domain has the same temperature and only the walls are heated, the simulation gets immediatelly unstable.

At the end of the simulation (with whatever tricks like bounding terms to get there) the simulation gets stable

Quote:

Originally Posted by agustinvo

The problems I found out with this model is the numericla stability: the leftSide as @jherb refers to can have a very bad condition number, which leads to spurious THF values, even if the region is isothermal. I'm looking for some workarounds for this.

Yuan Bao · November 27, 2021, 06:51

Quote:

Originally Posted by Fabio1893

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!

hi Fabio. I have the same question as you mentioned before. I have no idea to implement thermal expansion coefficient. I try to create a new value in createFields.H but failed.

zhenFoam · December 30, 2021, 02:08

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 jherb,

I couldn't find buoyantBoussinesqPimpleFoam in openfoam7. Maybe you mean openfoam6? I noticed that this FOAM has been deleted in openfoam7/8/9, I don't know why.

I also want to implant the AHFM-NRG model into OpenFOAM, just beginning, and I hope to get more guidance from you.

February 15, 2021, 14:03		#21
agustinvo Senior Member Agustín Villa Join Date: Apr 2013 Location: Alcorcón Posts: 313 Rep Power: 15	Hi @Andrea, nice to read you around! I managed to implement this some months ago in OpenFOAM 7, and it provides reasonable results in forced convection. As they say before, it seems quite unstable, so it is better to initialize with another turbulence model. If you are trying to use it in natural convection, it can be painful. About the fvSolutions, I use to relax the THF (check any pEqn.H to check out how you can do this), and try to use upwind schemes, at least at the beginning of the simulation. The problems I found out with this model is the numericla stability: the leftSide as @jherb refers to can have a very bad condition number, which leads to spurious THF values, even if the region is isothermal. I'm looking for some workarounds for this. Good luck! BTW: what is the test case you use to validate your implementation? Andrea1984 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 10:00
Wrong flow in ratating domain problem	Sanyo	CFX	17	August 15, 2015 06:20
CFX radiatin model wall heat flux imblance	cscfx	CFX	0	May 21, 2014 03:39
UDF for Heat Exchanger model	francois louw	FLUENT	2	July 16, 2010 02:21
Natural convection - Inlet boundary condition	max91	CFX	1	July 29, 2008 20:28

February 15, 2021, 14:15		#22
Andrea1984 Senior Member Andrea Join Date: Feb 2012 Location: Leeds, UK Posts: 179 Rep Power: 16	If it can be of any consolation, this turbulent heat flux closure was pretty unstable in both STAR-CCM+ and Code_Saturne as well. On top of what Augustin mentioned above, I would also suggest to ramp-up the values of the coefficients gradually from ~ 0 to their final value. It might also be helpful to set Ct4 = 0 altogether (unless the term multiplied by Ct4 plays a key role in the physics of the case that you are modelling). Good luck, Andrea