CFD Online Discussion Forums - Adding a source term in k-epsilon model

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:
https://image.ibb.co/f30FH0/1234.png
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!