k-epsilon implementation? (including volScalarField G)

sheaker · August 21, 2016, 12:24

Dear All.
I'm using: engineFoam (xiFoam), openfoam 2.2.1.

I'm trying to understand k-epsilon turbulence model in openFoam but I can't get through those k and epsilon equations.
Let's focus on k equation:

Code:

    tmp<fvScalarMatrix> kEqn
    (
        fvm::ddt(rho_, k_)
      + fvm::div(phi_, k_)
      - fvm::laplacian(DkEff(), k_)
     ==
        G
      - fvm::SuSp((2.0/3.0)*rho_*divU, k_)
      - fvm::Sp(rho_*epsilon_/k_, k_)
    );

I'm not sure if kEpsilon implements this equation but it looks similar:
$\frac{\partial\rho k}{\partial t}+\frac{\partial}{\partial x_i}(\rho V_i k) = \frac{\partial}{\partial x_j}[(\mu+\frac{\mu_{turb}}{\sigma_k})\frac{\partial k}{\partial x_j}]+\mu_{turb} (\sqrt{2S_{ij}S_{ij}})^2-\rho\epsilon$

If that's true, then:
fvm::ddt(rho_, k_) = $\frac{\partial\rho k}{\partial t}$
fvm::div(phi_, k_) = $\frac{\partial}{\partial x_i}(\rho V_i k)$
fvm::laplacian(DkEff(), k_) = $\frac{\partial}{\partial x_j}[(\mu+\frac{\mu_{turb}}{\sigma_k})\frac{\partial k}{\partial x_j}]$
fvm::Sp(rho_*epsilon_/k_, k_) = $\rho\epsilon$

Probably:
- fvm::SuSp((2.0/3.0)*rho_*divU, k_) = $-\frac{2}{3}*\rho*k*\frac{\partial V_i}{\partial x_i}$
and maybe this term is to remove trace but as long as I cant understand G, I can't go any further.

Definition of G is

Code:

 volScalarField G(GName(), mut_*(tgradU() && dev(twoSymm(tgradU()))));

and I don't know how to understand this.
I can't figure out how:
G - fvm::SuSp((2.0/3.0)*rho_*divU, k_) is equal to $\mu_{turb} (\sqrt{2S_{ij}S_{ij}})^2$

If someone can correct me or lead me to correct answer it would be great.
Greetings to all.

dradenkovic · August 23, 2016, 15:42

Have a look here (page 123, at the top)
http://www.tfd.chalmers.se/~lada/pos...-modelling.pdf

About G: You need volScalarField to enter k equation.

Regards,
Darko

sheaker · August 24, 2016, 16:17

Hello. Thank You for this pdf but I can't find full solution (exact implementation formula).

According to my further research k-epsilon turbulent model is based on this report:
"k-epsilon equation for compressible reciprocating engine flows", Journal of Energy, Vol. 7, No. 4 (1983), pp. 345-353.

Unfortunately it is not free.

dradenkovic · August 25, 2016, 05:34

Hello.

Here is your paper.

https://www.dropbox.com/s/iwnvw5om0b...48086.pdf?dl=0

About your first post, in exact k equation (have look here)

https://en.wikipedia.org/wiki/Turbulence_kinetic_energy

find production term. Insert Bussinesk assumption
$\overline{u_i^{\prime}u_j^{\prime}} = \frac{2}{3}\,k\delta_{ij} + ...$
in that production term and you will get terms that you need. In the case of incompressible flow,

$-\frac{2}{3}\,\rho\, k\,\frac{\partial V_i}{\partial x_i}$

is equal to zero.

Best regards,
Darko

sheaker · August 25, 2016, 14:21

Thank You once again, dear Darko. That 2nd report is very useful for me (eg. C_3 constant).

I'm really sorry but it still isn't clear for me.

I have done what you said.
The production term:

$P = -\rho \overline{u'_i u'_j}\frac{\partial \overline{u_i}}{\partial x_j}$

After applying Boussinesq Eddy Viscosity:
$P = - \frac{2}{3} \rho k \delta_{i,j}\frac{\partial \overline{u_i}}{\partial x_j}+\mu_{turb}(\frac{\partial \overline{u_i}}{\partial x_j}+\frac{\partial \overline{u_j}}{\partial x_i}) \frac{\partial \overline{u_i}}{\partial x_j}$

Removing delta Kronecker
$P = - \frac{2}{3} \rho k \frac{\partial \overline{u_i}}{\partial x_i}+\mu_{turb}(\frac{\partial \overline{u_i}}{\partial x_j}+\frac{\partial \overline{u_j}}{\partial x_i}) \frac{\partial \overline{u_i}}{\partial x_j}$

First term in openfoam:
- fvm::SuSp((2.0/3.0)*rho_*divU, k_) = $- \frac{2}{3} \rho k \frac{\partial \overline{u_i}}{\partial x_i}$

But I still don't understand how

G = mut_*(tgradU() && dev(twoSymm(tgradU()))) = $\mu_{turb}(\frac{\partial \overline{u_i}}{\partial x_j}+\frac{\partial \overline{u_j}}{\partial x_i}) \frac{\partial \overline{u_i}}{\partial x_j}$

tgradU is temporary field of gradient U

tgradU = $\frac{\partial \overline{u_i}}{\partial x_j}$
twoSymm(tgradU) returns 2*symmetric part of a tensor so:

twoSymm(tgradU) = $2* \frac{1}{2}(tgradU+tgradU^T) =$ $\frac{\partial \overline{u_i}}{\partial x_j}+\frac{\partial \overline{u_j}}{\partial x_i}$ (Am I right here? Im pretty sure but...)

dev(twoSymm(tgradU())) - this one returns only deviatoric component of tensor tgradU

dev(twoSymm(tgradU())) = $\frac{\partial \overline{u_i}}{\partial x_j}+\frac{\partial \overline{u_j}}{\partial x_i} - \frac{1}{3}(\frac{\partial \overline{u_k}}{\partial x_k}+\frac{\partial \overline{u_k}}{\partial x_k})\delta_{i,j} \mathbb{I}$
dev(twoSymm(tgradU())) = $\frac{\partial \overline{u_i}}{\partial x_j}+\frac{\partial \overline{u_j}}{\partial x_i} - \frac{2}{3}\frac{\partial \overline{u_k}}{\partial x_k}\delta_{i,j} \mathbb{I}$ Is that right?

And at last:
tgradU() && dev(twoSymm(tgradU())) means "double inner product" and it hard to understand.

Definition of double inner product for two 2nd rank tensors:
$A:B = A_{i,j}B_{k,l}\delta_{j,k} \delta_{i,l} = A_{i,j}B_{j,i}$

tgradU() && dev(twoSymm(tgradU())) = $[\frac{\partial \overline{u_i}}{\partial x_j}] : [\frac{\partial \overline{u_l}}{\partial x_k}+\frac{\partial \overline{u_k}}{\partial x_l} - \frac{2}{3}\frac{\partial \overline{u_m}}{\partial x_m}\delta_{k,l} \mathbb{I}]$

And I have no idea what to do now.

dradenkovic · August 26, 2016, 09:12

Book at this link has very good introduction part about tensors

https://www.dropbox.com/s/72rz5760b3...ou%29.pdf?dl=0

Sooner or later, you will have to learn it, so it is better to start now.
Pay attention to eq. 1.76, it will make things easier for you many times.

If we use (1.148) from above book, it follows:

tgradU() && dev(twoSymm(tgradU())) =
$\mu_t\,{\left({\mathbf{D} + \mathbf{\Omega}}\right)}:\left[2\mathbf{D} - \frac{2}{3}(tr \mathbf{D})\mathbf{I})\right] = \mu_t\left(2\mathbf{D}:\mathbf{D} - \frac{2}{3} (tr{\mathbf{D})^2}\right]$

Last term is equal to zero when flow is incompressible. It has been utilized that product of symmetric and asymmetric tensor is zero.

Regards,
Darko

sheaker · August 26, 2016, 18:26

Dear Darko. Thank You for a lot of patience to me.

I forgot about splitting gradU into symmetric and antisymmetric parts. I understand that:
(tgradU() && dev(twoSymm(tgradU())) = ${\left({\mathbf{D} + \mathbf{\Omega}}\right)}:\left[2\mathbf{D} - \frac{2}{3}(tr \mathbf{D})\mathbf{I}\right] = \left[2\mathbf{D}:\mathbf{D} - \frac{2}{3} (tr{\mathbf{D})^2}\right]$

According to previous hints:
$\left[{2\mathbf{D}:\mathbf{D} - \frac{2}{3} \left({tr{\mathbf{D}}}\right)^2}\right]$ should be equal to $\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)\frac{\partial u_i}{\partial x_j}$

And I'm not sure about that.
The first one is a scalar and it looks correct but I can't transform it to that second form.

Darko, wish You best!

dradenkovic · August 27, 2016, 06:25

I wasn't clear enough.

$\mathbf{D}$ is rate of strain tensor and from the above mentioned link from Chalmers, page 123, you can see that
${2\mathbf{D}}:\mathbf{D} =2 D_{ij}D_{ij} = \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right) \frac{\partial u_i}{\partial x_j}$

Term $- \frac{2}{3} \left(tr{\mathbf{D}}\right)^2$ is not in your starting equation.

When I think better from the very beginning we used assumption of incompressible flow through Bussinesq assumption - in the case of compressible flow, Bussinesq assumption is
$\mathbf{T} = 2\mu_t\left(\mathbf{D} - \frac{1}{3}(tr{\mathbf{D}}) \mathbf{I}\right) - \frac{2}{3} \rho k \mathbf{I}$

If you again insert Bussinesq assumption into production term of TKE, this forgotten term explains existence of $- \frac{2}{3} \left(tr{\mathbf{D}}\right)^2$ in Openfoam implementation of k-epsilon.

This agrees well with k equation of
http://turbmodels.larc.nasa.gov/ke-chien.html

Of course, in incompressible flow trace of rate of strain tensor is zero.

Regards,
Darko

sheaker · August 27, 2016, 10:07

Thank You for Your help. I think I get it now. So the therm $- \frac{2}{3} \left(tr{\mathbf{D}}\right)^2 = - \frac{2}{3} \left(\frac{\partial u_i}{\partial x_i}\right)^2$ SHOULD exist in my TKE equation?

$\frac{\partial\rho k}{\partial t}+\frac{\partial}{\partial x_i}(\rho u_i k) = \frac{\partial}{\partial x_j}$ $\left[ \left( \mu+\frac{\mu_{turb}}{\sigma_k} \right) \frac{\partial k}{\partial x_j} \right]+\mu_{turb} \left[ \left( \frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i} \right) \frac{\partial u_i}{\partial x_j} - \frac{2}{3} \left(\frac{\partial u_i}{\partial x_i}\right)^2 \right] - \frac{2}{3}\rho k\frac{\partial u_i}{\partial x_i}-\rho\epsilon$

And G(openFoam) = $\mu_{turb} \left[ \left( \frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i} \right) \frac{\partial u_i}{\partial x_j} - \frac{2}{3} \left(\frac{\partial u_i}{\partial x_i}\right)^2 \right]$
Production P = $\mu_{turb} \left[ \left( \frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i} \right) \frac{\partial u_i}{\partial x_j} - \frac{2}{3} \left(\frac{\partial u_i}{\partial x_i}\right)^2 \right] - \frac{2}{3}\rho k\frac{\partial u_i}{\partial x_i}$

I think that is correct now.

Kindest regards,
Oskar

dradenkovic · August 27, 2016, 12:46

Couple hours ago I verified this on link in my above post (but now, at the moment of this writing, above link doesn't work).

If you check eq. (24) and eq. (37) from paper "k-e Equation for Compressible Reciprocating Engine Flows", 1983, I believe you will see that is it.

Regards,
Darko

Mehrez · May 22, 2018, 10:51

Hello all,
Thank you for this very useful discussion.
I'm using multiphaseEulerFoam (OFv4) to perform simulation of two-phase (gas dispersed in liquid) dispersed flow. To model turbulence, I used the k-epsilon model.
The code has a little bit changed in version 4 (https://github.com/OpenFOAM/OpenFOAM...lon/kEpsilon.C) compared to what you shared in last posts.
In OFv4, we have source terms "kSource()" and "epsilonSource()". My understanding is that they are used for bubble induced turbulence in the case of two-phase dispersed flow.
Do you know how they are implemented and what is their exact definition?
Thank you.
mhrz

calf.Z · March 22, 2019, 08:02

I am also interested in kSource（）and epsilonSource（). What is its function and meaning?

jherb · July 1, 2019, 19:29

They are used to modify the turbulence model, e.g. for buoyancy effects. See https://github.com/OpenFOAM/OpenFOAM...uoyantKEpsilon

This turbulence model inherits from the basic compressible k-epsilon turbulence model and overwrites those methods to implement the effect of buoyancy on turbulence.

Quote:

Originally Posted by calf.Z

I am also interested in kSource（）and epsilonSource（). What is its function and meaning?

sourav8016 · May 20, 2021, 10:51

Quote:

Originally Posted by dradenkovic

Hello.

Here is your paper.

https://www.dropbox.com/s/iwnvw5om0b...48086.pdf?dl=0

About your first post, in exact k equation (have look here)

https://en.wikipedia.org/wiki/Turbulence_kinetic_energy

find production term. Insert Bussinesk assumption
$\overline{u_i^{\prime}u_j^{\prime}} = \frac{2}{3}\,k\delta_{ij} + ...$
in that production term and you will get terms that you need. In the case of incompressible flow,

$-\frac{2}{3}\,\rho\, k\,\frac{\partial V_i}{\partial x_i}$

is equal to zero.

Best regards,
Darko

Dear Sir, could you please upload the paper again? The link is showing that the paper is deleted.

sourav8016 · May 20, 2021, 11:03

Quote:

Originally Posted by sheaker

Thank You for Your help. I think I get it now. So the therm $- \frac{2}{3} \left(tr{\mathbf{D}}\right)^2 = - \frac{2}{3} \left(\frac{\partial u_i}{\partial x_i}\right)^2$ SHOULD exist in my TKE equation?

$\frac{\partial\rho k}{\partial t}+\frac{\partial}{\partial x_i}(\rho u_i k) = \frac{\partial}{\partial x_j}$ $\left[ \left( \mu+\frac{\mu_{turb}}{\sigma_k} \right) \frac{\partial k}{\partial x_j} \right]+\mu_{turb} \left[ \left( \frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i} \right) \frac{\partial u_i}{\partial x_j} - \frac{2}{3} \left(\frac{\partial u_i}{\partial x_i}\right)^2 \right] - \frac{2}{3}\rho k\frac{\partial u_i}{\partial x_i}-\rho\epsilon$

And G(openFoam) = $\mu_{turb} \left[ \left( \frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i} \right) \frac{\partial u_i}{\partial x_j} - \frac{2}{3} \left(\frac{\partial u_i}{\partial x_i}\right)^2 \right]$
Production P = $\mu_{turb} \left[ \left( \frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i} \right) \frac{\partial u_i}{\partial x_j} - \frac{2}{3} \left(\frac{\partial u_i}{\partial x_i}\right)^2 \right] - \frac{2}{3}\rho k\frac{\partial u_i}{\partial x_i}$

I think that is correct now.

Kindest regards,
Oskar

Dear sir, Could you please share the paper and the book shared by Darko Radenkovi? The link is broken.

kooki_13 · April 6, 2022, 06:54

Dear Sheaker

Do you have any Idea about G?

Here we have
tmp<fvScalarMatrix> kEqn
(
fvm::ddt(alpha, rho, k_)
+ fvm::div(alphaRhoPhi, k_)
- fvm::laplacian(alpha*rho*DkEff(), k_)
==
alpha*rho*G
- fvm::SuSp((2.0/3.0)*alpha*rho*divU, k_)
- fvm::Sp(Ce(D, KK)*alpha*rho*sqrt(k_)/this->delta(), k_)
+ kSource()
+ fvModels.source(alpha, rho, k_)
);

I have understood the term - fvm::SuSp((2.0/3.0)*alpha*rho*divU, k_)
but still this one is unclear:
alpha*rho*G

kooki_13 · April 6, 2022, 06:56

Quote:

Originally Posted by sheaker

Thank You for Your help. I think I get it now. So the therm $- \frac{2}{3} \left(tr{\mathbf{D}}\right)^2 = - \frac{2}{3} \left(\frac{\partial u_i}{\partial x_i}\right)^2$ SHOULD exist in my TKE equation?

$\frac{\partial\rho k}{\partial t}+\frac{\partial}{\partial x_i}(\rho u_i k) = \frac{\partial}{\partial x_j}$ $\left[ \left( \mu+\frac{\mu_{turb}}{\sigma_k} \right) \frac{\partial k}{\partial x_j} \right]+\mu_{turb} \left[ \left( \frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i} \right) \frac{\partial u_i}{\partial x_j} - \frac{2}{3} \left(\frac{\partial u_i}{\partial x_i}\right)^2 \right] - \frac{2}{3}\rho k\frac{\partial u_i}{\partial x_i}-\rho\epsilon$

And G(openFoam) = $\mu_{turb} \left[ \left( \frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i} \right) \frac{\partial u_i}{\partial x_j} - \frac{2}{3} \left(\frac{\partial u_i}{\partial x_i}\right)^2 \right]$
Production P = $\mu_{turb} \left[ \left( \frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i} \right) \frac{\partial u_i}{\partial x_j} - \frac{2}{3} \left(\frac{\partial u_i}{\partial x_i}\right)^2 \right] - \frac{2}{3}\rho k\frac{\partial u_i}{\partial x_i}$

I think that is correct now.

Kindest regards,
Oskar

Dear Sheaker

Do you have any Idea about G?

Here we have
tmp<fvScalarMatrix> kEqn
(
fvm::ddt(alpha, rho, k_)
+ fvm::div(alphaRhoPhi, k_)
- fvm::laplacian(alpha*rho*DkEff(), k_)
==
alpha*rho*G
- fvm::SuSp((2.0/3.0)*alpha*rho*divU, k_)
- fvm::Sp(Ce(D, KK)*alpha*rho*sqrt(k_)/this->delta(), k_)
+ kSource()
+ fvModels.source(alpha, rho, k_)
);

I have understood the term - fvm::SuSp((2.0/3.0)*alpha*rho*divU, k_)
but still this one is unclear:
alpha*rho*G

mechkween · May 31, 2022, 11:27

G refers to the production of the TKE. This is written out as:

Code:

nut*(tgradU() && dev(twoSymm(tgradU())))

which is essentially a double inner over gradU(j,i) and 2*S(i,j) - 2/3 S(k,k) \delta(i,j). The first term is the velocity gradient and the second term represents the model for the Reynolds stress less the nut, which is pre-multiplied.

alainislas · June 6, 2023, 20:41

Dear Sourav Hossain

The book is "Papanastasiou, T., Georgiou, G., & Alexandrou, A. N. (2021). Viscous fluid flow. CRC press."

August 21, 2016, 12:24	k-epsilon implementation? (including volScalarField G)	#1
sheaker Senior Member Oskar Join Date: Nov 2015 Location: Poland Posts: 184 Rep Power: 11	Dear All. I'm using: engineFoam (xiFoam), openfoam 2.2.1. I'm trying to understand k-epsilon turbulence model in openFoam but I can't get through those k and epsilon equations. Let's focus on k equation: Code: tmp<fvScalarMatrix> kEqn ( fvm::ddt(rho_, k_) + fvm::div(phi_, k_) - fvm::laplacian(DkEff(), k_) == G - fvm::SuSp((2.0/3.0)rho_divU, k_) - fvm::Sp(rho_epsilon_/k_, k_) ); I'm not sure if kEpsilon implements this equation but it looks similar: $\frac{\partial\rho k}{\partial t}+\frac{\partial}{\partial x_i}(\rho V_i k) = \frac{\partial}{\partial x_j}[(\mu+\frac{\mu_{turb}}{\sigma_k})\frac{\partial k}{\partial x_j}]+\mu_{turb} (\sqrt{2S_{ij}S_{ij}})^2-\rho\epsilon$ If that's true, then: fvm::ddt(rho_, k_) = $\frac{\partial\rho k}{\partial t}$ fvm::div(phi_, k_) = $\frac{\partial}{\partial x_i}(\rho V_i k)$ fvm::laplacian(DkEff(), k_) = $\frac{\partial}{\partial x_j}[(\mu+\frac{\mu_{turb}}{\sigma_k})\frac{\partial k}{\partial x_j}]$ fvm::Sp(rho_epsilon_/k_, k_) = $\rho\epsilon$ Probably: - fvm::SuSp((2.0/3.0)rho_divU, k_) = $-\frac{2}{3}\rhok\frac{\partial V_i}{\partial x_i}$ and maybe this term is to remove trace but as long as I cant understand G, I can't go any further. Definition of G is Code: volScalarField G(GName(), mut_(tgradU() && dev(twoSymm(tgradU())))); and I don't know how to understand this. I can't figure out how: G - fvm::SuSp((2.0/3.0)rho_divU, k_) is equal to $\mu_{turb} (\sqrt{2S_{ij}S_{ij}})^2$ If someone can correct me or lead me to correct answer it would be great. Greetings to all. randolph, Cajal and pramodhkumaresh like this.

August 23, 2016, 15:42		#2
dradenkovic Member Darko Radenkovic Join Date: Oct 2015 Posts: 38 Rep Power: 11	Have a look here (page 123, at the top) http://www.tfd.chalmers.se/~lada/pos...-modelling.pdf About G: You need volScalarField to enter k equation. Regards, Darko sheaker likes this.

August 25, 2016, 05:34		#4
dradenkovic Member Darko Radenkovic Join Date: Oct 2015 Posts: 38 Rep Power: 11	Hello. Here is your paper. https://www.dropbox.com/s/iwnvw5om0b...48086.pdf?dl=0 About your first post, in exact k equation (have look here) https://en.wikipedia.org/wiki/Turbulence_kinetic_energy find production term. Insert Bussinesk assumption $\overline{u_i^{\prime}u_j^{\prime}} = \frac{2}{3}\,k\delta_{ij} + ...$ in that production term and you will get terms that you need. In the case of incompressible flow, $-\frac{2}{3}\,\rho\, k\,\frac{\partial V_i}{\partial x_i}$ is equal to zero. Best regards, Darko sheaker likes this.

August 25, 2016, 14:21		#5
sheaker Senior Member Oskar Join Date: Nov 2015 Location: Poland Posts: 184 Rep Power: 11	Thank You once again, dear Darko. That 2nd report is very useful for me (eg. C_3 constant). I'm really sorry but it still isn't clear for me. I have done what you said. The production term: $P = -\rho \overline{u'_i u'_j}\frac{\partial \overline{u_i}}{\partial x_j}$ After applying Boussinesq Eddy Viscosity: $P = - \frac{2}{3} \rho k \delta_{i,j}\frac{\partial \overline{u_i}}{\partial x_j}+\mu_{turb}(\frac{\partial \overline{u_i}}{\partial x_j}+\frac{\partial \overline{u_j}}{\partial x_i}) \frac{\partial \overline{u_i}}{\partial x_j}$ Removing delta Kronecker $P = - \frac{2}{3} \rho k \frac{\partial \overline{u_i}}{\partial x_i}+\mu_{turb}(\frac{\partial \overline{u_i}}{\partial x_j}+\frac{\partial \overline{u_j}}{\partial x_i}) \frac{\partial \overline{u_i}}{\partial x_j}$ First term in openfoam: - fvm::SuSp((2.0/3.0)rho_divU, k_) = $- \frac{2}{3} \rho k \frac{\partial \overline{u_i}}{\partial x_i}$ But I still don't understand how G = mut_(tgradU() && dev(twoSymm(tgradU()))) = $\mu_{turb}(\frac{\partial \overline{u_i}}{\partial x_j}+\frac{\partial \overline{u_j}}{\partial x_i}) \frac{\partial \overline{u_i}}{\partial x_j}$ tgradU is temporary field of gradient U tgradU = $\frac{\partial \overline{u_i}}{\partial x_j}$ twoSymm(tgradU) returns 2symmetric part of a tensor so: twoSymm(tgradU) = $2* \frac{1}{2}(tgradU+tgradU^T) =$ $\frac{\partial \overline{u_i}}{\partial x_j}+\frac{\partial \overline{u_j}}{\partial x_i}$ (Am I right here? Im pretty sure but...) dev(twoSymm(tgradU())) - this one returns only deviatoric component of tensor tgradU dev(twoSymm(tgradU())) = $\frac{\partial \overline{u_i}}{\partial x_j}+\frac{\partial \overline{u_j}}{\partial x_i} - \frac{1}{3}(\frac{\partial \overline{u_k}}{\partial x_k}+\frac{\partial \overline{u_k}}{\partial x_k})\delta_{i,j} \mathbb{I}$ dev(twoSymm(tgradU())) = $\frac{\partial \overline{u_i}}{\partial x_j}+\frac{\partial \overline{u_j}}{\partial x_i} - \frac{2}{3}\frac{\partial \overline{u_k}}{\partial x_k}\delta_{i,j} \mathbb{I}$ Is that right? And at last: tgradU() && dev(twoSymm(tgradU())) means "double inner product" and it hard to understand. Definition of double inner product for two 2nd rank tensors: $A:B = A_{i,j}B_{k,l}\delta_{j,k} \delta_{i,l} = A_{i,j}B_{j,i}$ tgradU() && dev(twoSymm(tgradU())) = $[\frac{\partial \overline{u_i}}{\partial x_j}] : [\frac{\partial \overline{u_l}}{\partial x_k}+\frac{\partial \overline{u_k}}{\partial x_l} - \frac{2}{3}\frac{\partial \overline{u_m}}{\partial x_m}\delta_{k,l} \mathbb{I}]$ And I have no idea what to do now. randolph likes this. Last edited by sheaker; August 25, 2016 at 18:18.

August 26, 2016, 09:12		#6
dradenkovic Member Darko Radenkovic Join Date: Oct 2015 Posts: 38 Rep Power: 11	Book at this link has very good introduction part about tensors https://www.dropbox.com/s/72rz5760b3...ou%29.pdf?dl=0 Sooner or later, you will have to learn it, so it is better to start now. Pay attention to eq. 1.76, it will make things easier for you many times. If we use (1.148) from above book, it follows: tgradU() && dev(twoSymm(tgradU())) = $\mu_t\,{\left({\mathbf{D} + \mathbf{\Omega}}\right)}:\left[2\mathbf{D} - \frac{2}{3}(tr \mathbf{D})\mathbf{I})\right] = \mu_t\left(2\mathbf{D}:\mathbf{D} - \frac{2}{3} (tr{\mathbf{D})^2}\right]$ Last term is equal to zero when flow is incompressible. It has been utilized that product of symmetric and asymmetric tensor is zero. Regards, Darko randolph and sheaker like this.

August 24, 2016, 16:17		#3
sheaker Senior Member Oskar Join Date: Nov 2015 Location: Poland Posts: 184 Rep Power: 11	Hello. Thank You for this pdf but I can't find full solution (exact implementation formula). According to my further research k-epsilon turbulent model is based on this report: "k-epsilon equation for compressible reciprocating engine flows", Journal of Energy, Vol. 7, No. 4 (1983), pp. 345-353. Unfortunately it is not free.

August 26, 2016, 18:26		#7
sheaker Senior Member Oskar Join Date: Nov 2015 Location: Poland Posts: 184 Rep Power: 11	Dear Darko. Thank You for a lot of patience to me. I forgot about splitting gradU into symmetric and antisymmetric parts. I understand that: (tgradU() && dev(twoSymm(tgradU())) = ${\left({\mathbf{D} + \mathbf{\Omega}}\right)}:\left[2\mathbf{D} - \frac{2}{3}(tr \mathbf{D})\mathbf{I}\right] = \left[2\mathbf{D}:\mathbf{D} - \frac{2}{3} (tr{\mathbf{D})^2}\right]$ According to previous hints: $\left[{2\mathbf{D}:\mathbf{D} - \frac{2}{3} \left({tr{\mathbf{D}}}\right)^2}\right]$ should be equal to $\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)\frac{\partial u_i}{\partial x_j}$ And I'm not sure about that. The first one is a scalar and it looks correct but I can't transform it to that second form. Darko, wish You best! Last edited by sheaker; August 27, 2016 at 05:32.

August 27, 2016, 06:25		#8
dradenkovic Member Darko Radenkovic Join Date: Oct 2015 Posts: 38 Rep Power: 11	I wasn't clear enough. $\mathbf{D}$ is rate of strain tensor and from the above mentioned link from Chalmers, page 123, you can see that ${2\mathbf{D}}:\mathbf{D} =2 D_{ij}D_{ij} = \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right) \frac{\partial u_i}{\partial x_j}$ Term $- \frac{2}{3} \left(tr{\mathbf{D}}\right)^2$ is not in your starting equation. When I think better from the very beginning we used assumption of incompressible flow through Bussinesq assumption - in the case of compressible flow, Bussinesq assumption is $\mathbf{T} = 2\mu_t\left(\mathbf{D} - \frac{1}{3}(tr{\mathbf{D}}) \mathbf{I}\right) - \frac{2}{3} \rho k \mathbf{I}$ If you again insert Bussinesq assumption into production term of TKE, this forgotten term explains existence of $- \frac{2}{3} \left(tr{\mathbf{D}}\right)^2$ in Openfoam implementation of k-epsilon. This agrees well with k equation of http://turbmodels.larc.nasa.gov/ke-chien.html Of course, in incompressible flow trace of rate of strain tensor is zero. Regards, Darko

August 27, 2016, 10:07		#9
sheaker Senior Member Oskar Join Date: Nov 2015 Location: Poland Posts: 184 Rep Power: 11	Thank You for Your help. I think I get it now. So the therm $- \frac{2}{3} \left(tr{\mathbf{D}}\right)^2 = - \frac{2}{3} \left(\frac{\partial u_i}{\partial x_i}\right)^2$ SHOULD exist in my TKE equation? $\frac{\partial\rho k}{\partial t}+\frac{\partial}{\partial x_i}(\rho u_i k) = \frac{\partial}{\partial x_j}$ $\left[ \left( \mu+\frac{\mu_{turb}}{\sigma_k} \right) \frac{\partial k}{\partial x_j} \right]+\mu_{turb} \left[ \left( \frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i} \right) \frac{\partial u_i}{\partial x_j} - \frac{2}{3} \left(\frac{\partial u_i}{\partial x_i}\right)^2 \right] - \frac{2}{3}\rho k\frac{\partial u_i}{\partial x_i}-\rho\epsilon$ And G(openFoam) = $\mu_{turb} \left[ \left( \frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i} \right) \frac{\partial u_i}{\partial x_j} - \frac{2}{3} \left(\frac{\partial u_i}{\partial x_i}\right)^2 \right]$ Production P = $\mu_{turb} \left[ \left( \frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i} \right) \frac{\partial u_i}{\partial x_j} - \frac{2}{3} \left(\frac{\partial u_i}{\partial x_i}\right)^2 \right] - \frac{2}{3}\rho k\frac{\partial u_i}{\partial x_i}$ I think that is correct now. Kindest regards, Oskar randolph likes this.

August 27, 2016, 12:46		#10
dradenkovic Member Darko Radenkovic Join Date: Oct 2015 Posts: 38 Rep Power: 11	Couple hours ago I verified this on link in my above post (but now, at the moment of this writing, above link doesn't work). If you check eq. (24) and eq. (37) from paper "k-e Equation for Compressible Reciprocating Engine Flows", 1983, I believe you will see that is it. Regards, Darko sheaker likes this.

May 22, 2018, 10:51		#11
Mehrez Member Sami Join Date: Nov 2012 Location: Cap Town, South Africa Posts: 87 Rep Power: 14	Hello all, Thank you for this very useful discussion. I'm using multiphaseEulerFoam (OFv4) to perform simulation of two-phase (gas dispersed in liquid) dispersed flow. To model turbulence, I used the k-epsilon model. The code has a little bit changed in version 4 (https://github.com/OpenFOAM/OpenFOAM...lon/kEpsilon.C) compared to what you shared in last posts. In OFv4, we have source terms "kSource()" and "epsilonSource()". My understanding is that they are used for bubble induced turbulence in the case of two-phase dispersed flow. Do you know how they are implemented and what is their exact definition? Thank you. mhrz

March 22, 2019, 08:02		#12
calf.Z Senior Member Jianrui Zeng Join Date: May 2018 Location: China Posts: 157 Rep Power: 8	I am also interested in kSource（）and epsilonSource（). What is its function and meaning?

April 6, 2022, 06:54		#16
kooki_13 New Member Join Date: Jul 2017 Posts: 14 Rep Power: 9	Dear Sheaker Do you have any Idea about G? Here we have tmp<fvScalarMatrix> kEqn ( fvm::ddt(alpha, rho, k_) + fvm::div(alphaRhoPhi, k_) - fvm::laplacian(alpharhoDkEff(), k_) == alpharhoG - fvm::SuSp((2.0/3.0)alpharhodivU, k_) - fvm::Sp(Ce(D, KK)alpharhosqrt(k_)/this->delta(), k_) + kSource() + fvModels.source(alpha, rho, k_) ); I have understood the term - fvm::SuSp((2.0/3.0)alpharhodivU, k_) but still this one is unclear: alpharho*G

May 31, 2022, 11:27		#18
mechkween Member Join Date: Dec 2012 Posts: 33 Rep Power: 14	G refers to the production of the TKE. This is written out as: Code: nut(tgradU() && dev(twoSymm(tgradU()))) which is essentially a double inner over gradU(j,i) and 2S(i,j) - 2/3 S(k,k) \delta(i,j). The first term is the velocity gradient and the second term represents the model for the Reynolds stress less the nut, which is pre-multiplied. reverseila likes this.

June 6, 2023, 20:41		#19
alainislas Senior Member Alain Islas Join Date: Nov 2019 Location: Mexico Posts: 142 Rep Power: 7	Dear Sourav Hossain The book is "Papanastasiou, T., Georgiou, G., & Alexandrou, A. N. (2021). Viscous fluid flow. CRC press."

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
SimpleFoam k and epsilon bounded	nedved	OpenFOAM Running, Solving & CFD	16	March 4, 2017 09:30
[mesh manipulation] Importing Multiple Meshes	thomasnwalshiii	OpenFOAM Meshing & Mesh Conversion	18	December 19, 2015 19:57
Creating a new field from terms of the turbulence model	HaZe	OpenFOAM Programming & Development	15	November 24, 2014 14:51
writing execFlowFunctionObjects	immortality	OpenFOAM Post-Processing	30	September 15, 2013 07:16
SimpleFoam k and epsilon bounded	nedved	OpenFOAM Running, Solving & CFD	1	November 25, 2008 21:21