/*********************************************************************/
/*Secondary Settling Tank Modeling								*/
/*Solids Concentration UDS										*/
/*Date: 22/10/05 										    		*/
/*Authors: 
	Roberto Garcia (roberto.garcia@fluent.fr)
	Didier Bessette (didier.bessette@fluent.fr)
	Benoit Post (benoit.post@fluent.fr)    							*/
/*********************************************************************/
	
#include "udf.h"
#include "turb.h"

/* Schmidt number */
#define SC_T 0.7

/* Material properties (SI units) */
#define RHOP 1450. /* dry particle density */
#define RHOW 1000. /* clean water density */

/* Gravity (SI units) */
#define GX 9.81 /* gravity force (m/s-2) */

/* Tackacs settling velocity */
/* US = US0 x exp[-RH(C-Cmin)] - US0 x exp[-RP(C-Cmin)] */
/* Constants (SI units) */
#define US0 0.005	/* reference settling velocity (m/s) */
#define RH 0.7		/*floc settling parameter (m3/kg)*/
#define RP 5		/* colloidal settling parameter (m3/kg) */
#define CMIN 0.01	/* nonsettleable concentration (kg/m3) */
#define USMAX 0.002	/* maximum settling velocity (m/s) */

static real settling_velocity(real C)
{
  real US;

  if(C<CMIN)
    US=0.;
  else
    US=US0*(exp(-RH*(C-CMIN))-exp(-RP*(C-CMIN)));

  if(US>0.002)
    US=0.002;

  return US;
}

/* Convective UDS concentration term */
DEFINE_UDS_FLUX(settling_flux,f,tf,i)
{
  cell_t c0=F_C0(f,tf),c1=F_C1(f,tf);
  Thread *t0=THREAD_T0(tf),*t1=THREAD_T1(tf);
  real face_concentration;
  real gravity_vector[]={1,0,0},area[ND_ND];/* Gravity vector (1,0,0)*/

  if ((t1==NULL)||(THREAD_TYPE(tf)==THREAD_F_WALL))
    /* if this is a boundary zone or a 0-thickness wall */
    face_concentration=F_UDSI(f,tf,0);
  else
    /* if this an interior face, the face value is the average of the upstream and downstream cells */
    face_concentration=.5*(C_UDSI(c0,t0,0)+C_UDSI(c1,t1,0));
  /* F_AREA macro stores the components (Ax,Ay) of the area vector on the array area */
  F_AREA(area,f,tf);

  if (THREAD_TYPE(tf)==THREAD_F_WALL)
    return 0.;
  else
    return F_FLUX(f,tf)+C_R(c0,t0)*settling_velocity(face_concentration)*NV_DOT(gravity_vector,area);
	/*if the face is parallel to the gravity_vector, NV_DOT(gravity_vector,area) = 0 */
	/*if the face is perpendicular to the gravity_vector, NV_DOT(gravity_vector,area) = 1 */
}

/* Diffusivity UDS concentration term */
DEFINE_DIFFUSIVITY(turbulent_diff,c,tc,i)
{
  /* UDS Diffusivity is defined as the turbulent dynamic viscoisty over the Schmidt number */
  return 1e-10+C_MU_T(c,tc)/SC_T;
}

/* Gravity source term to account for particle-fluid density differences */
DEFINE_SOURCE(X_mom_src,c,tc,ds,eqn)
{
  /* UDM-0 retrieves settling velocity and it is accesible during standard post-processing */
  C_UDMI(c,tc,0)=settling_velocity(C_UDSI(c,tc,0));

  /* UDM-1 retrieves UDS Diffusivity term */  
  C_UDMI(c,tc,1)=C_MU_T(c,tc)/SC_T;
  
  /* UDM-2 retrieves gravity source term and it is accesible during standard post-processing */
  C_UDMI(c,tc,2) = GX*C_UDSI(c,tc,0)*(RHOP-C_R(c,tc))/RHOP;
  
  return C_UDMI(c,tc,2) ;
}

/* Turbulent kinetic energy equation source term */
DEFINE_SOURCE(turb_k_source,c,tc,ds,eqn)
{
  real beta=(RHOP-RHOW)/(RHOP*RHOW);
  /* UDM-3 retrieves k source term and it is accesible during standard post-processing */  
  /* C_UDSI_G accesses the UDS gradient */
  C_UDMI(c,tc,3)=-GX*beta*C_MU_T(c,tc)/SC_T*C_UDSI_G(c,tc,0)[0];
  ds[eqn]=-2*GX*beta*0.09*C_K(c,tc)/(C_D(c,tc)*SC_T)*C_UDSI_G(c,tc,0)[0];
  
  return C_UDMI(c,tc,3);
}

/* Turbulent energy dissipation equation source term */
DEFINE_SOURCE(turb_e_source,c,tc,ds,eqn)
{
  real C3_eps=0.;
  /* C3_eps is calculated according to the equation 11.4-26 of Fluent manual */
  C3_eps=tanh(fabs(C_U(c,tc)/C_V(c,tc))+fabs(SMALL_VP));
  /* UDM-4 retrieves eps source term and it is accesible during standard post-processing */
  C_UDMI(c,tc,4)=C3_eps*1.44*C_D(c,tc)/C_K(c,tc)*C_UDMI(c,tc,3);

  return C_UDMI(c,tc,4);
}

/* Viscosity is defined as specified in reference: 								*/
/* Computing Shear Flow and Sludge Blanket in Secondary Clarifiers 				*/
/* Djamel Lakehal, Peter Krebs, Johan Krijgsman, and Wolfgang Rodi, Member, ASCE	*/
DEFINE_PROPERTY(viscosity,c,tc)
{
  if (C_UDSI(c,tc,0)>0.7)
    return 0.00327*pow(10.,0.132*C_UDSI(c,tc,0));
  else
    if (C_UDSI(c,tc,0)>0.01)
      return 0.00404527222+0.00304527222*(C_UDSI(c,tc,0)-0.7)/0.69;
  else
    return 0.001;
}
