Adding a source term in k-epsilon model

cshsgy · November 17, 2018, 03: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, 07:22

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

November 17, 2018, 03:14	Adding a source term in k-epsilon model	#1
cshsgy New Member Chen Sihe Join Date: Sep 2017 Posts: 26 Rep Power: 10	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!

Thread Tools	Search this Thread
Show Printable Version Email this Page	Search this Thread: Advanced Search
Display Modes
Linear Mode Switch to Hybrid Mode Switch to Threaded Mode

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
[swak4Foam] swak4foam for OpenFOAM 4.0	mnikku	OpenFOAM Community Contributions	80	May 17, 2022 09: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 19:13
Problem compiling a custom Lagrangian library	brbbhatti	OpenFOAM Programming & Development	2	July 7, 2014 12:32
OpenFOAM without MPI	kokizzu	OpenFOAM Installation	4	May 26, 2014 10:17
Trouble compiling utilities using source-built OpenFOAM	Artur	OpenFOAM Programming & Development	14	October 29, 2013 11:59

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