Adding a source term in k-epsilon model

cshsgy · November 17, 2018, 02:14

Hi all,
I am testing out a model raised by S.B. Pope, which modified the epsilon equation with a source term. I tried to implement it with udf in the FLUENT, but the convergence seems to be a problem. The modified epsilon-equation looks like the following:
${\frac {\partial (\rho \epsilon )}{\partial t}}+{\frac {\partial (\rho \epsilon u_{i})}{\partial x_{i}}}={\frac {\partial }{\partial x_{j}}}\left[{\frac {\mu _{t}}{\sigma _{\epsilon }}}{\frac {\partial \epsilon }{\partial x_{j}}}\right]+C_{1\epsilon }{\frac {\epsilon }{k}}2{\mu _{t}}{E_{ij}}{E_{ij}}-(C_{2\epsilon }-C_{3\epsilon}\chi)\rho {\frac {\epsilon ^{2}}{k}}$
The added source term is related to the $C_{3\epsilon}\chi$ term, which is a measure of the vortex stretching. Here, $\chi=\frac{1}{4}(\frac{k}{\epsilon})^3(\frac{\partial U}{\partial r}-\frac{\partial V}{\partial x})^2\frac{V}{r}$ .
The case I am simulating is basically a turbulent round jet. I firstly calculated the case until convergence, and then implemented the source term to the epsilon equation as the following code:

Code:

#include "udf.h"
DEFINE_SOURCE(epsilon_Pope_source, cell, thread, dS, eqn)
{
real source;
real time_scale = C_K(cell,thread) / C_D(cell,thread);
real omega = 0.5 * time_scale * (C_DUDY(cell,thread)-C_DVDX(cell,thread));
real x_c[ND_ND];
real kai;
C_CENTROID(x_c,cell,thread);
kai = time_scale * omega * omega * C_V(cell,thread)/x_c[1];// Use DVDY instead of V/r?

source = 0.79 * kai / time_scale * C_D(cell,thread);
dS[eqn] = source * (-1) / C_D(cell,thread);
return source;
}

However, the convergence is not possible. The plot of residual is as follows:

Clearly the white line (continuity) is driven to a very high level. Also other variables would not converge. Interestingly if I do not initiate with the calculation of the standard model, the residual stabilizes at a even higher value. Cannot figure out what is going wrong.
I tried to observe the flow field, and find out that the v,k,epsilon fields are all making no sense. I think the problem comes with the udf and I might have not done it properly. What I am doing is simply to add the $C_{\epsilon3}\chi\frac{\epsilon^2}{k}$ term as the source term directly. What is the issue with that?
Thanks so much for your help!

andresLearns · November 2, 2021, 06:22

Hi cshsgy, did you solve the issue? I plan to do a similar implementation.

November 17, 2018, 02:14	Adding a source term in k-epsilon model	#1
cshsgy New Member Chen Sihe Join Date: Sep 2017 Posts: 26 Rep Power: 8	Hi all, I am testing out a model raised by S.B. Pope, which modified the epsilon equation with a source term. I tried to implement it with udf in the FLUENT, but the convergence seems to be a problem. The modified epsilon-equation looks like the following: ${\frac {\partial (\rho \epsilon )}{\partial t}}+{\frac {\partial (\rho \epsilon u_{i})}{\partial x_{i}}}={\frac {\partial }{\partial x_{j}}}\left[{\frac {\mu _{t}}{\sigma _{\epsilon }}}{\frac {\partial \epsilon }{\partial x_{j}}}\right]+C_{1\epsilon }{\frac {\epsilon }{k}}2{\mu _{t}}{E_{ij}}{E_{ij}}-(C_{2\epsilon }-C_{3\epsilon}\chi)\rho {\frac {\epsilon ^{2}}{k}}$ The added source term is related to the $C_{3\epsilon}\chi$ term, which is a measure of the vortex stretching. Here, $\chi=\frac{1}{4}(\frac{k}{\epsilon})^3(\frac{\partial U}{\partial r}-\frac{\partial V}{\partial x})^2\frac{V}{r}$ . The case I am simulating is basically a turbulent round jet. I firstly calculated the case until convergence, and then implemented the source term to the epsilon equation as the following code: Code: #include "udf.h" DEFINE_SOURCE(epsilon_Pope_source, cell, thread, dS, eqn) { real source; real time_scale = C_K(cell,thread) / C_D(cell,thread); real omega = 0.5 * time_scale * (C_DUDY(cell,thread)-C_DVDX(cell,thread)); real x_c[ND_ND]; real kai; C_CENTROID(x_c,cell,thread); kai = time_scale * omega * omega * C_V(cell,thread)/x_c[1];// Use DVDY instead of V/r? source = 0.79 * kai / time_scale * C_D(cell,thread); dS[eqn] = source * (-1) / C_D(cell,thread); return source; } However, the convergence is not possible. The plot of residual is as follows: Clearly the white line (continuity) is driven to a very high level. Also other variables would not converge. Interestingly if I do not initiate with the calculation of the standard model, the residual stabilizes at a even higher value. Cannot figure out what is going wrong. I tried to observe the flow field, and find out that the v,k,epsilon fields are all making no sense. I think the problem comes with the udf and I might have not done it properly. What I am doing is simply to add the $C_{\epsilon3}\chi\frac{\epsilon^2}{k}$ term as the source term directly. What is the issue with that? Thanks so much for your help!

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
[swak4Foam] swak4foam for OpenFOAM 4.0	mnikku	OpenFOAM Community Contributions	80	May 17, 2022 08:06
[foam-extend.org] Problems installing foam-extend-4.0 on openSUSE 42.2 and Ubuntu 16.04	ordinary	OpenFOAM Installation	19	September 3, 2019 18:13
Problem compiling a custom Lagrangian library	brbbhatti	OpenFOAM Programming & Development	2	July 7, 2014 11:32
OpenFOAM without MPI	kokizzu	OpenFOAM Installation	4	May 26, 2014 09:17
Trouble compiling utilities using source-built OpenFOAM	Artur	OpenFOAM Programming & Development	14	October 29, 2013 10:59

November 2, 2021, 06:22		#2
andresLearns New Member Andres Join Date: Aug 2021 Posts: 3 Rep Power: 4	Hi cshsgy, did you solve the issue? I plan to do a similar implementation.