/* turbulence model: v2f by Durbin (1991) and modified by Kalizin (1999)				      */
/* UDF (User Defined Function) Routine for FLUENT                                                             */
/* see also the paper:                                                                                        */
/* Implementation and comparison of different turbulence models for three dimensional wall jets with Fluent   */
/* FLUENT CFD Forum 2005 Bad Nauheim, Germany								      */
/* Heschl Ch*., Sanz W.**, Klanatsky P.*                                                                      */
/* *  Fachhochschulstudiengänge Burgenland GmbH, Steinamangerstraße 21, A-7423 Pinkafeld                      */
/* ** TU Graz, Institut of Thermal Turbomachinery and Machine Dynamics, Inffeldgasse 25, A-8010 Graz          */
/* Contact: christian.heschl@fh-burgenland.at, wolfgang.sanz@tugraz.at, peter.klanatsky@fh-burgenland.at      */

/* Die Autoren übernhemen keine Verantwortung oder Haftung für den Inhalt der unten angeführten UDF-Routinen  */


#include "udf.h"
#include "math.h"

/* Turbulence model constants */
const real C_MU=0.22;
const real CE_D=0.050;
const real CE_1=1.44;
const real CE_2=1.92;
const real SIG_K=1.0;
const real SIG_E=1.3;
const real C_1=1.4;
const real C_2=0.3;
const real C_L=0.23;
const real C_ETA=70.0;
const real y_star_limit=30.0;
const real C_MU_ke=0.09;

/* User-defined scalars */
enum
{
K,
E,
V2,
F,
N_REQUIRED_UDS
};




DEFINE_SOURCE(k_source, c, t, dS, eqn)
{
real source;


dS[eqn]=- C_R(c,t)/C_UDMI(c,t,1);
source=C_UDMI(c,t,3) - C_R(c,t)*C_UDSI(c,t,K)/C_UDMI(c,t,1);
return source;
}



DEFINE_SOURCE(e_source, c, t, dS, eqn)
{
//real CE_1;
real source;


//CE_1=1.4 * (1.0+CE_D*sqrt( C_UDSI(c,t,K) / C_UDSI(c,t,V2) ));

dS[eqn]= -CE_2*C_R(c,t)/C_UDMI(c,t,1);
source= CE_1*C_UDMI(c,t,3)/C_UDMI(c,t,1)  -  CE_2*C_R(c,t)*C_UDSI(c,t,E)/C_UDMI(c,t,1);
return source;
}



DEFINE_SOURCE(v2_source, c, t, dS, eqn)
{
real source;
real kf;


kf=C_UDSI(c,t,K)*C_UDSI(c,t,F);

dS[eqn]=-6.*C_R(c,t)*C_UDSI(c,t,E)/C_UDSI(c,t,K);
source=C_R(c,t)*kf-6.*C_R(c,t)*C_UDSI(c,t,V2)*C_UDSI(c,t,E)/C_UDSI(c,t,K);
return source;
}



DEFINE_SOURCE(f_source, c, t, dS, eqn)
{
real source;


dS[eqn]=-1.0/SQR(C_UDMI(c,t,2));
source=-C_UDSI(c,t,F)/SQR(C_UDMI(c,t,2))-1.0/(SQR(C_UDMI(c,t,2))*C_UDMI(c,t,1))*((C_1-6.)*C_UDSI(c,t,V2)/C_UDSI(c,t,K)-2./3.*(C_1-1.))+C_2*C_UDMI(c,t,3)/(C_R(c,t)*C_UDSI(c,t,K)*SQR(C_UDMI(c,t,2)));
return source;
}

DEFINE_DIFFUSIVITY(ke_v2f_diffusivity, c, t, eqn)
{
real diff;

switch (eqn)
{
case K:
diff=C_UDMI(c,t,0)/SIG_K+C_MU_L(c,t);
break;

case E:
diff=C_UDMI(c,t,0)/SIG_E+C_MU_L(c,t);
break;

case V2:
diff=C_UDMI(c,t,0)+C_MU_L(c,t);
break;

case F:
diff=1.0;
break;

default:
diff=0.0;
}
return diff;
}

DEFINE_UDS_FLUX(user_flux, f, t, eqn)
{

switch (eqn)
{
case K:
return F_FLUX(f,t);
break;

case E:
return F_FLUX(f,t);
break;

case V2:
return F_FLUX(f,t);
break;

case F:
return 0.0;
break;

default:
return 0.0;
}

}



DEFINE_ADJUST(ke_adjust,domain)
{

	Thread *t;
	cell_t c;
	real y_star, L, T, V;
	real mu_t;

	/* Set the turbulent viscosity */
	thread_loop_c(t,domain)
		if (FLUID_THREAD_P(t))
		{
			begin_c_loop(c,t)
			{				

			y_star=C_R(c,t)*pow(C_MU_ke,0.25)*sqrt(C_UDSI(c,t,K))*C_WALL_DIST(c,t)/C_MU_L(c,t);

			if (y_star <= y_star_limit)
			{
				T=MAX( C_UDSI(c,t,K)/C_UDSI(c,t,E), 6.0*sqrt( C_MU_L(c,t)/(C_R(c,t)*C_UDSI(c,t,E))) );
				L=C_L*MAX(pow( (C_UDSI(c,t,K)),(3./2.))/C_UDSI(c,t,E),C_ETA*pow( ((pow(C_MU_L(c,t)/C_R(c,t),3.)/C_UDSI(c,t,E))),0.25));
			}
			else
			{
				T=C_UDSI(c,t,K)/C_UDSI(c,t,E);
				L=C_L*pow( C_UDSI(c,t,K),(3./2.))/C_UDSI(c,t,E);
			}
	
				mu_t=C_R(c,t)*C_MU*C_UDSI(c,t,V2)*T;

				C_K(c,t)=C_UDSI(c,t,K);
				C_D(c,t)=C_UDSI(c,t,E);

				C_UDMI(c,t,0)=mu_t;
				C_UDMI(c,t,1)=T;
				C_UDMI(c,t,2)=L;
				C_UDMI(c,t,3)=mu_t*SQR(Strainrate_Mag(c,t));

			}
			end_c_loop(c,t)
		}
}



DEFINE_TURBULENT_VISCOSITY(user_mu_t,c,t)
{

return C_UDMI(c,t,0);
}



DEFINE_PROFILE(e_bc, t, position)
{
real dy;

face_t f;
cell_t c0;
Thread *t0=t->t0;
real xw[ND_ND], xc[ND_ND], dx[ND_ND];

begin_f_loop(f,t)
{
c0=F_C0(f,t);
F_CENTROID(xw,f,t);
C_CENTROID(xc,c0,t0);
NV_VV(dx, =, xc, -, xw);
dy=ND_MAG(dx[0], dx[1], dx[2]);

F_PROFILE(f,t,position)=2.*C_MU_L(c0,t0)/C_R(c0,t0)*C_UDSI(c0,t0,K)/SQR(dy);

}
end_f_loop(f,t)
}