CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

Implementation of turbulent algebraic heat flux model

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

Like Tree3Likes

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   January 8, 2021, 06:34
Default Implementation of turbulent algebraic heat flux model
  #1
New Member
 
Join Date: Jan 2021
Posts: 10
Rep Power: 5
Fabio1893 is on a distinguished road
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!
Minghao_Li likes this.

Last edited by Fabio1893; January 8, 2021 at 07:38.
Fabio1893 is offline   Reply With Quote

Old   January 11, 2021, 08:51
Default
  #2
Senior Member
 
Joachim Herb
Join Date: Sep 2010
Posts: 650
Rep Power: 21
jherb is on a distinguished road
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 is offline   Reply With Quote

Old   January 11, 2021, 11:02
Default
  #3
New Member
 
Join Date: Jan 2021
Posts: 10
Rep Power: 5
Fabio1893 is on a distinguished road
Quote:
Originally Posted by jherb View Post
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!

Last edited by Fabio1893; January 12, 2021 at 08:54.
Fabio1893 is offline   Reply With Quote

Old   January 12, 2021, 08:59
Default
  #4
New Member
 
Join Date: Jan 2021
Posts: 10
Rep Power: 5
Fabio1893 is on a distinguished road
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?

Last edited by Fabio1893; January 12, 2021 at 10:05.
Fabio1893 is offline   Reply With Quote

Old   January 16, 2021, 13:08
Default
  #5
Senior Member
 
Joachim Herb
Join Date: Sep 2010
Posts: 650
Rep Power: 21
jherb is on a distinguished road
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.
Farid likes this.
jherb is offline   Reply With Quote

Old   January 16, 2021, 13:25
Default
  #6
Senior Member
 
Joachim Herb
Join Date: Sep 2010
Posts: 650
Rep Power: 21
jherb is on a distinguished road
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
jherb is offline   Reply With Quote

Old   January 19, 2021, 04:37
Default
  #7
New Member
 
Join Date: Jan 2021
Posts: 10
Rep Power: 5
Fabio1893 is on a distinguished road
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.
Fabio1893 is offline   Reply With Quote

Old   January 19, 2021, 15:11
Default
  #8
Senior Member
 
Joachim Herb
Join Date: Sep 2010
Posts: 650
Rep Power: 21
jherb is on a distinguished road
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).
jherb is offline   Reply With Quote

Old   January 20, 2021, 03:42
Default
  #9
New Member
 
Join Date: Jan 2021
Posts: 10
Rep Power: 5
Fabio1893 is on a distinguished road
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?
Fabio1893 is offline   Reply With Quote

Old   January 22, 2021, 04:13
Default
  #10
Senior Member
 
Andrea
Join Date: Feb 2012
Location: Leeds, UK
Posts: 179
Rep Power: 16
Andrea1984 is on a distinguished road
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 View Post
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?
Andrea1984 is offline   Reply With Quote

Old   January 22, 2021, 15:38
Default
  #11
Senior Member
 
Joachim Herb
Join Date: Sep 2010
Posts: 650
Rep Power: 21
jherb is on a distinguished road
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.
jherb is offline   Reply With Quote

Old   January 23, 2021, 11:20
Default
  #12
Senior Member
 
Andrea
Join Date: Feb 2012
Location: Leeds, UK
Posts: 179
Rep Power: 16
Andrea1984 is on a distinguished road
No problem, I fully understand your point.

BTW I have been working at NRG on this stuff for a while
Andrea1984 is offline   Reply With Quote

Old   January 26, 2021, 06:23
Default
  #13
New Member
 
Join Date: Jan 2021
Posts: 10
Rep Power: 5
Fabio1893 is on a distinguished road
Quote:
Originally Posted by jherb View Post
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
Fabio1893 is offline   Reply With Quote

Old   January 29, 2021, 17:50
Default
  #14
Senior Member
 
Joachim Herb
Join Date: Sep 2010
Posts: 650
Rep Power: 21
jherb is on a distinguished road
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 View Post
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 is offline   Reply With Quote

Old   February 2, 2021, 08:53
Default
  #15
New Member
 
Join Date: Jan 2021
Posts: 10
Rep Power: 5
Fabio1893 is on a distinguished road
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
Fabio1893 is offline   Reply With Quote

Old   February 2, 2021, 18:07
Default
  #16
Senior Member
 
Joachim Herb
Join Date: Sep 2010
Posts: 650
Rep Power: 21
jherb is on a distinguished road
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 is offline   Reply With Quote

Old   February 2, 2021, 18:28
Default
  #17
Senior Member
 
Joachim Herb
Join Date: Sep 2010
Posts: 650
Rep Power: 21
jherb is on a distinguished road
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.
jherb is offline   Reply With Quote

Old   February 4, 2021, 06:06
Default
  #18
New Member
 
Join Date: Jan 2021
Posts: 10
Rep Power: 5
Fabio1893 is on a distinguished road
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
Fabio1893 is offline   Reply With Quote

Old   February 7, 2021, 11:04
Default
  #19
Senior Member
 
Joachim Herb
Join Date: Sep 2010
Posts: 650
Rep Power: 21
jherb is on a distinguished road
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.
jherb is offline   Reply With Quote

Old   February 11, 2021, 09:24
Default
  #20
New Member
 
Join Date: Jan 2021
Posts: 10
Rep Power: 5
Fabio1893 is on a distinguished road
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
Fabio1893 is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


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


All times are GMT -4. The time now is 08:58.