/* turbulence model: nonlinear v2f by Petersson Reif (1999), v2f 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;
const real zeta=8.0/3.0;

/* 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=MIN( C_UDSI(c,t,K)*C_UDSI(c,t,F) , -1./C_UDMI(c,t,1)*((C_1-6.)*C_UDSI(c,t,V2)-2./3.*C_UDSI(c,t,K)*(C_1-1.))+C_2*C_UDMI(c,t,3)/C_R(c,t) );

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, T1;
	real mu_t_ke, mu_t_v2f, mu_t;
	real S11, S12, S13, S21, S22, S23, S31, S32, S33;
	real W11, W12, W13, W21, W22, W23, W31, W32, W33;
	real S, W, eta1, eta2, eta3;
	real beta1, gamma1, cmu1, cmu2, cmu3;

	/* 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_v2f= C_R(c,t)*C_MU*C_UDSI(c,t,V2)*T;
				mu_t_ke = C_R(c,t)*C_MU_ke*SQR(C_UDSI(c,t,K))/C_UDSI(c,t,E);
				mu_t=MIN(mu_t_v2f,mu_t_ke);

				V=MAX( 2./3.-C_UDSI(c,t,V2)/C_UDSI(c,t,K) , 0.);

				S11 = 0.5*( C_DUDX(c,t)+C_DUDX(c,t) );
				S12 = 0.5*( C_DUDY(c,t)+C_DVDX(c,t) );
				S13 = 0.5*( C_DUDZ(c,t)+C_DWDX(c,t) );
				S21 = 0.5*( C_DVDX(c,t)+C_DUDY(c,t) );
				S22 = 0.5*( C_DVDY(c,t)+C_DVDY(c,t) );
				S23 = 0.5*( C_DVDZ(c,t)+C_DWDY(c,t) );
				S31 = 0.5*( C_DWDX(c,t)+C_DUDZ(c,t) );
				S32 = 0.5*( C_DWDY(c,t)+C_DVDZ(c,t) );
				S33 = 0.5*( C_DWDZ(c,t)+C_DWDZ(c,t) );
				W11 = 0.5*( C_DUDX(c,t)-C_DUDX(c,t) );
				W12 = 0.5*( C_DUDY(c,t)-C_DVDX(c,t) );
				W13 = 0.5*( C_DUDZ(c,t)-C_DWDX(c,t) );
				W21 = 0.5*( C_DVDX(c,t)-C_DUDY(c,t) );
				W22 = 0.5*( C_DVDY(c,t)-C_DVDY(c,t) );
				W23 = 0.5*( C_DVDZ(c,t)-C_DWDY(c,t) );
				W31 = 0.5*( C_DWDX(c,t)-C_DUDZ(c,t) );
				W32 = 0.5*( C_DWDY(c,t)-C_DVDZ(c,t) );
				W33 = 0.5*( C_DWDZ(c,t)-C_DWDZ(c,t) );

				//S=S11*S11+S12*S21+S13*S31+S21*S12+S22*S22+S23*S32+S31*S13+S32*S23+S33*S33;
				//W=-(W11*W11+W12*W21+W13*W31+W21*W12+W22*W22+W23*W32+W31*W13+W32*W23+W33*W33);
				S=S11*S11+S12*S12+S13*S13+S21*S21+S22*S22+S23*S23+S31*S31+S32*S32+S33*S33;
				W=W11*W11+W12*W12+W13*W13+W21*W21+W22*W22+W23*W23+W31*W31+W32*W32+W33*W33;
	
				eta1=SQR(T)*S;
				eta2=SQR(T)*W;
				eta3=eta1-eta2;

				beta1=1./(0.1+sqrt(eta1*eta2));
				gamma1=1./(0.1+eta1);

				cmu1=C_MU;
				cmu2=6.0/5.0*sqrt(MAX(1.0-SQR(cmu1*C_UDSI(c,t,V2)/C_UDSI(c,t,K))*2.*eta1,0.))/(beta1+sqrt(eta1*eta2));
				cmu3=6.0/5.0/(gamma1+eta1);

				T1=MIN( T , 2./3.*C_UDSI(c,t,K)/(cmu1*C_UDSI(c,t,V2)*sqrt(S*zeta)) );

								
				C_RUU(c,t)=	-2.0*cmu1*C_UDSI(c,t,V2)*T1*S11 + 2./3.*C_UDSI(c,t,K)
							-V*C_UDSI(c,t,K)*SQR(T)*
							( cmu2*(S11*W11+S12*W21+S13*W31+S11*W11+S12*W21+S13*W31) - cmu3*(S11*S11+S12*S12+S13*S13-1./3.*S) );
				
				C_RVV(c,t)=-2.0*cmu1*C_UDSI(c,t,V2)*T1*S22 + 2./3.*C_UDSI(c,t,K)
							-V*C_UDSI(c,t,K)*SQR(T)*
							( cmu2*(S21*W12+S22*W22+S23*W32+S21*W12+S22*W22+S23*W32) - cmu3*(S12*S12+S22*S22+S23*S23-1./3.*S) );


				C_RWW(c,t)= -2.0*cmu1*C_UDSI(c,t,V2)*T1*S33 + 2./3.*C_UDSI(c,t,K)
							-V*C_UDSI(c,t,K)*SQR(T)*
							( cmu2*(S31*W13+S32*W23+S33*W33+S31*W13+S32*W23+S33*W33) - cmu3*(S13*S13+S23*S23+S33*S33-1./3.*S) );
				
				C_RUV(c,t)= -2.0*cmu1*C_UDSI(c,t,V2)*T1*S12
							-V*C_UDSI(c,t,K)*SQR(T)*
							( cmu2*(S11*W12+S12*W22+S13*W32+S21*W11+S22*W21+S23*W31) - cmu3*(S11*S12+S12*S22+S13*S32) );

				C_RUW(c,t)= -2.0*cmu1*C_UDSI(c,t,V2)*T1*S13
							-V*C_UDSI(c,t,K)*SQR(T)*
							( cmu2*(S11*W13+S12*W23+S13*W33+S31*W11+S32*W21+S33*W31) - cmu3*(S11*S13+S12*S23+S13*S33) );

				C_RVW(c,t)= -2.0*cmu1*C_UDSI(c,t,V2)*T1*S23
							-V*C_UDSI(c,t,K)*SQR(T)*
							( cmu2*(S21*W13+S22*W23+S23*W33+S31*W12+S32*W22+S33*W32) - cmu3*(S21*S13+S22*S23+S23*S33) );


				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_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)
}


DEFINE_ON_DEMAND(on_demand_calc)
{
Domain *d; /* declare domain pointer since it is not passed as an argument to the DEFINE macro */

Thread *t;
cell_t c;

d = Get_Domain(1); /* Get the domain using Fluent utility */

thread_loop_c(t,d)
{
begin_c_loop(c,t)
{
				C_UDSI(c,t,K) = C_K(c,t);
				C_UDSI(c,t,E) = C_D(c,t);
				C_UDSI(c,t,V2) = 2.0/3.0*C_K(c,t);
				C_UDSI(c,t,F)=0.01;
}
end_c_loop(c,t)
}
}