wallfunction.c File Reference

#include "wallfunction.h"
#include <math.h>
#include <stdlib.h>
#include <stdio.h>

Include dependency graph for wallfunction.c:

Go to the source code of this file.

Macros
#define	KAPPA   0.41
#define	LOGLAW_B   5.5

Functions
void	noslip (UserCtx *user, double sc, double sb, Cmpnts Ua, Cmpnts Uc, Cmpnts *Ub, double nx, double ny, double nz)
void	freeslip (UserCtx *user, double sc, double sb, Cmpnts Ua, Cmpnts Uc, Cmpnts *Ub, double nx, double ny, double nz)
double	E_coeff (double utau, double ks, double nu)
double	u_hydset_roughness (double nu, double y, double utau, double ks)
double	f_hydset (double nu, double u, double y, double utau0, double ks)
double	df_hydset (double nu, double u, double y, double utau0, double ks)
double	find_utau_hydset (double nu, double u, double y, double utau_guess, double ks)
double	nu_t (double yplus)
double	integrate_1 (double nu, double y, double utau, int m)
double	taw (double nu, double utau, double y, double u, double dpdt)
double	u_Cabot (double nu, double y, double utau, double dpdt, double taw)
double	u_Werner (double nu, double y, double utau)
double	f_Werner (double nu, double u, double y, double utau)
double	df_Werner (double nu, double u, double y, double utau)
double	f_Cabot (double nu, double u, double y, double utau, double dpdt, double dpdtn)
double	df_Cabot (double nu, double u, double y, double utau, double dpdt, double dpdtn)
void	find_utau_Cabot (double nu, double u, double y, double guess, double dpdt, double dpdtn, double *utau, double *taw1, double *taw2)
double	find_utau_Werner (double nu, double u, double y, double guess)
double	sign (double a)
double	u_loglaw (double y, double utau, double roughness)
double	find_utau_loglaw (double u, double y, double roughness)
void	wall_function (UserCtx *user, double sc, double sb, Cmpnts Ua, Cmpnts Uc, Cmpnts *Ub, PetscReal *ustar, double nx, double ny, double nz)
void	wall_function_loglaw (UserCtx *user, double ks, double sc, double sb, Cmpnts Ua, Cmpnts Uc, Cmpnts *Ub, PetscReal *ustar, double nx, double ny, double nz)
void	wall_function_Cabot (UserCtx *user, double ks, double sc, double sb, Cmpnts Ua, Cmpnts Uc, Cmpnts *Ub, PetscReal *ustar, double nx, double ny, double nz, double dpdx, double dpdy, double dpdz, int count)

Macro Definition Documentation

KAPPA

#define KAPPA   0.41

Definition at line 10 of file wallfunction.c.

LOGLAW_B

#define LOGLAW_B   5.5

Definition at line 11 of file wallfunction.c.

Function Documentation

noslip()

void noslip	(	UserCtx *	user,
		double	sc,
		double	sb,
		Cmpnts	Ua,
		Cmpnts	Uc,
		Cmpnts *	Ub,
		double	nx,
		double	ny,
		double	nz
	)

Definition at line 13 of file wallfunction.c.

{
    double u_c = Uc.x - Ua.x, v_c = Uc.y - Ua.y, w_c = Uc.z - Ua.z;
    double un = u_c * nx + v_c * ny + w_c * nz;
    double ut = u_c - un * nx, vt = v_c - un * ny, wt = w_c - un * nz;
    double ut_mag = sqrt( ut*ut + vt*vt + wt*wt );
    
    (*Ub).x = sb/sc * u_c;
    (*Ub).y = sb/sc * v_c;
    (*Ub).z = sb/sc * w_c;
    
    (*Ub).x += Ua.x;
    (*Ub).y += Ua.y;
    (*Ub).z += Ua.z;
}

Here is the caller graph for this function:

freeslip()

void freeslip	(	UserCtx *	user,
		double	sc,
		double	sb,
		Cmpnts	Ua,
		Cmpnts	Uc,
		Cmpnts *	Ub,
		double	nx,
		double	ny,
		double	nz
	)

Definition at line 30 of file wallfunction.c.

{
  //  in the Normal direction interpolation 
  double Uan = Ua.x * nx + Ua.y * ny + Ua.z * nz;
  double Ucn = Uc.x * nx + Uc.y * ny + Uc.z * nz;
  double Ubn = Uan + (Ucn - Uan) * sb/sc;
  //   in the Tangential direction Ub_t = Uc_t
  double ut = Uc.x - Ucn * nx, vt = Uc.y - Ucn * ny, wt = Uc.z - Ucn * nz;
  
    (*Ub).x = ut;
    (*Ub).y = vt;
    (*Ub).z = wt;
    
    (*Ub).x += Ubn*nx;
    (*Ub).y += Ubn*ny;
    (*Ub).z += Ubn*nz;
}

E_coeff()

double E_coeff	(	double	utau,
		double	ks,
		double	nu
	)

Definition at line 50 of file wallfunction.c.

{
  double kplus=utau*ks/nu, dB;
  if(kplus<=2.25) dB = 0.0;
  else if(kplus>2.25 && kplus<90.) dB = (LOGLAW_B-8.5+1./KAPPA*log(kplus))*sin(0.4258*(log(kplus)-0.811));
  else if(kplus>=90.) dB = LOGLAW_B-8.5+1./KAPPA*log(kplus);
  return exp(KAPPA*(LOGLAW_B-dB));
}

Here is the caller graph for this function:

u_hydset_roughness()

double u_hydset_roughness	(	double	nu,
		double	y,
		double	utau,
		double	ks
	)

Definition at line 59 of file wallfunction.c.

{
  double y0plus=11.53, y1plus=300,f;
  double yplus = utau*y/nu;
 
  if(yplus<=y0plus){f= utau * yplus;}
  else if (yplus<=y1plus) {f= utau/KAPPA*log(E_coeff(utau,ks,nu)*yplus);}
  else {f=-1.;}
  return f;
}

Here is the call graph for this function:

Here is the caller graph for this function:

f_hydset()

double f_hydset	(	double	nu,
		double	u,
		double	y,
		double	utau0,
		double	ks
	)

Definition at line 70 of file wallfunction.c.

{
double y0plus=11.53, f;
double yplus=utau0*y/nu;
 if (yplus<=y0plus) {f= utau0*yplus-u;}
 else {f= utau0*(1./KAPPA*log(E_coeff(utau0,ks,nu)*yplus))-u;}
return f;
}

Here is the call graph for this function:

Here is the caller graph for this function:

df_hydset()

double df_hydset	(	double	nu,
		double	u,
		double	y,
		double	utau0,
		double	ks
	)

Definition at line 79 of file wallfunction.c.

{
double eps=1.e-7;
return (f_hydset(nu, u, y, utau0 + eps, ks) - f_hydset(nu, u, y, utau0 - eps, ks))/(2.*eps);
}

Here is the call graph for this function:

Here is the caller graph for this function:

find_utau_hydset()

double find_utau_hydset	(	double	nu,
		double	u,
		double	y,
		double	utau_guess,
		double	ks
	)

Definition at line 86 of file wallfunction.c.

{
 double utau,utau0 = utau_guess;  
 int ii;
 for(ii=0;ii<30;ii++){
 utau=utau0-f_hydset(nu,u,y,utau0,ks)/df_hydset(nu,u,y,utau0,ks);
 if (fabs(utau-utau0)<1.e-7)break;
 utau0=utau;
  }
 if (ii==30) LOG_ALLOW(LOCAL,LOG_DEBUG, "iter utau diff %d %le %le  %le\n",ii,utau,utau0,fabs(utau-utau0));
 //printf("iter utau diff %d %le %le  %le\n",ii,utau,utau0,fabs(utau-utau0));
return utau;
};

Here is the call graph for this function:

Here is the caller graph for this function:

nu_t()

double nu_t

(

double

yplus

)

Definition at line 103 of file wallfunction.c.

{
    return KAPPA * yplus * pow ( 1. - exp( - yplus / 19. ) , 2.0 );
};

Here is the caller graph for this function:

integrate_1()

double integrate_1	(	double	nu,
		double	y,
		double	utau,
		int	m
	)

Definition at line 107 of file wallfunction.c.

{
    int ny=30;  //number of integartion points
    int i;
    double nut,f=0.0,yplus;
    double dy=y/(ny-1);
    double F[ny];
    for (i=0;i<ny;i++){
        yplus=(i*dy)*utau/nu;
        nut=nu_t(yplus);
        if (m==0){
         F[i]=1.0/(1.0+nut)/nu;
        }   else
        {
            F[i]=dy*i/(1.0+nut)/nu;
        }
    }
    for (i=1;i<ny-1;i++){
        f+=(F[i-1]+2*F[i]+F[i+1])*0.25*dy;
    }
    f+=.1667*dy*(2*F[0]+F[1]+2*F[ny-1]+F[ny-2]);
    return f;
}

Here is the call graph for this function:

Here is the caller graph for this function:

taw()

double taw	(	double	nu,
		double	utau,
		double	y,
		double	u,
		double	dpdt
	)

Definition at line 130 of file wallfunction.c.

{
    double yplus = y*utau/nu;
    double sum=0.0;
    double f1 = integrate_1(nu,y,utau,0);
    double f2 = integrate_1(nu,y,utau,1);
    return 1/f1*(u-dpdt*f2); 
}

Here is the call graph for this function:

Here is the caller graph for this function:

u_Cabot()

double u_Cabot	(	double	nu,
		double	y,
		double	utau,
		double	dpdt,
		double	taw
	)

Definition at line 138 of file wallfunction.c.

{
    double f1 = integrate_1(nu,y,utau,0);
    double f2 = integrate_1(nu,y,utau,1);
    return  taw*f1+dpdt*f2;
};

Here is the call graph for this function:

Here is the caller graph for this function:

u_Werner()

double u_Werner	(	double	nu,
		double	y,
		double	utau
	)

Definition at line 145 of file wallfunction.c.

{
    double yplus = utau * y / nu; 
    double A=8.3, B=1./7.;
    
    if(yplus<=11.81) return yplus*utau;
    else return A * pow(yplus, B) * utau;
};

Here is the caller graph for this function:

f_Werner()

double f_Werner	(	double	nu,
		double	u,
		double	y,
		double	utau
	)

Definition at line 154 of file wallfunction.c.

{
    double ym=11.81*nu/utau;
    double A=8.3, B=1./7.;
    
    if( fabs(u) <= nu/(2*ym) * pow(A, 2./(1.-B) ) ) return utau*utau - u/y*nu;
    else return utau*utau -  u/fabs(u) * pow( 0.5*(1-B) * pow(A, (1+B)/(1-B)) * pow(nu/y, 1+B) + (1+B)/A * pow(nu/y, B) * fabs(u), 2/(1+B) );
}

Here is the caller graph for this function:

df_Werner()

double df_Werner	(	double	nu,
		double	u,
		double	y,
		double	utau
	)

Definition at line 163 of file wallfunction.c.

{
    double eps=1.e-7;
    return ( f_Werner(nu, u, y, utau+eps) - f_Werner(nu, u, y, utau-eps) ) / ( 2*eps ) ;
}

Here is the call graph for this function:

Here is the caller graph for this function:

f_Cabot()

double f_Cabot	(	double	nu,
		double	u,
		double	y,
		double	utau,
		double	dpdt,
		double	dpdtn
	)

Definition at line 169 of file wallfunction.c.

{
    return utau - sqrt(fabs(taw( nu, utau, y, u, dpdt )));
}

Here is the call graph for this function:

Here is the caller graph for this function:

df_Cabot()

double df_Cabot	(	double	nu,
		double	u,
		double	y,
		double	utau,
		double	dpdt,
		double	dpdtn
	)

Definition at line 174 of file wallfunction.c.

{
    double eps=1.e-7;
    return ( f_Cabot(nu, u, y, utau+eps, dpdt,dpdtn) - f_Cabot(nu, u, y, utau-eps, dpdt, dpdtn)) / ( 2*eps ) ;
}

Here is the call graph for this function:

Here is the caller graph for this function:

find_utau_Cabot()

void find_utau_Cabot	(	double	nu,
		double	u,
		double	y,
		double	guess,
		double	dpdt,
		double	dpdtn,
		double *	utau,
		double *	taw1,
		double *	taw2
	)

Definition at line 191 of file wallfunction.c.

{
    double x, x0=guess;
    
    int i;
    
    for(i=0; i<30; i++) {
        x =fabs (x0 - f_Cabot(nu, u, y, x0, dpdt, dpdtn)/df_Cabot(nu, u, y, x0, dpdt, dpdtn));
        if( fabs(x0 - x) < 1.e-10) break;
        x0 = x;
    }
    
    if( fabs(x0 - x) > 1.e-5 && i>=29 ) printf("\n!!!!!!!! Iteration Failed !!!!!!!!\n");
        *utau=x;
        *taw1=taw( nu, x, y, u, dpdt );
        *taw2=taw( nu, x, y, 0.0, dpdtn );
 
};

Here is the call graph for this function:

Here is the caller graph for this function:

find_utau_Werner()

double find_utau_Werner	(	double	nu,
		double	u,
		double	y,
		double	guess
	)

Definition at line 210 of file wallfunction.c.

{
    double x, x0=guess;
    
    int i;
    
    for(i=0; i<20; i++) {
        x = x0 - f_Werner(nu, u, y, x0)/df_Werner(nu, u, y, x0);
        if( fabs(x0 - x) < 1.e-7 ) break;
        x0 = x;
    }
    
    if( fabs(x0 - x) > 1.e-5 && i>=19 ) printf("\n!!!!!!!! Iteration Failed !!!!!!!!\n");
        
    return x;
};

Here is the call graph for this function:

sign()

double sign

(

double

a

)

Definition at line 227 of file wallfunction.c.

{
    if(a>0) return 1;
    else if(a<0) return -1;
    else return 0;
}

u_loglaw()

double u_loglaw	(	double	y,
		double	utau,
		double	roughness
	)

Definition at line 235 of file wallfunction.c.

{
    return utau *  1./KAPPA * log( (roughness+y) / roughness ) ;
}

find_utau_loglaw()

double find_utau_loglaw	(	double	u,
		double	y,
		double	roughness
	)

Definition at line 240 of file wallfunction.c.

{
    return KAPPA * u / log( (y+roughness) / roughness);
};

wall_function()

void wall_function	(	UserCtx *	user,
		double	sc,
		double	sb,
		Cmpnts	Ua,
		Cmpnts	Uc,
		Cmpnts *	Ub,
		PetscReal *	ustar,
		double	nx,
		double	ny,
		double	nz
	)

Definition at line 245 of file wallfunction.c.

{
    SimCtx *simCtx = user->simCtx;
        double u_c = Uc.x - Ua.x, v_c = Uc.y - Ua.y, w_c = Uc.z - Ua.z;
    double un = u_c * nx + v_c * ny + w_c * nz;
    double ut = u_c - un * nx, vt = v_c - un * ny, wt = w_c - un * nz;
    double ut_mag = sqrt( ut*ut + vt*vt + wt*wt );
    //*ustar = find_utau_Cabot(1./simCtx->ren, ut_mag, sc, 0.01, 0);
    double utau=0.05;
    double ut_mag_modeled = u_Werner(1./simCtx->ren, sb, utau);
    //double ut_mag_modeled = u_Cabot(1./simCtx->ren, sb, *ustar, 0);
    
    if(ut_mag>1.e-10) {
        ut *= ut_mag_modeled/ut_mag;
        vt *= ut_mag_modeled/ut_mag;
        wt *= ut_mag_modeled/ut_mag;
    }
    else ut=vt=wt=0;
                    
    // u = ut + (u.n)n
    (*Ub).x = ut + sb/sc * un * nx;
    (*Ub).y = vt + sb/sc * un * ny;
    (*Ub).z = wt + sb/sc * un * nz;
    
    (*Ub).x += Ua.x;
    (*Ub).y += Ua.y;
    (*Ub).z += Ua.z;
}

Here is the call graph for this function:

Here is the caller graph for this function:

wall_function_loglaw()

void wall_function_loglaw	(	UserCtx *	user,
		double	ks,
		double	sc,
		double	sb,
		Cmpnts	Ua,
		Cmpnts	Uc,
		Cmpnts *	Ub,
		PetscReal *	ustar,
		double	nx,
		double	ny,
		double	nz
	)

Definition at line 276 of file wallfunction.c.

{
    SimCtx *simCtx = user->simCtx;
        double u_c = Uc.x - Ua.x, v_c = Uc.y - Ua.y, w_c = Uc.z - Ua.z;
    //  double ua_n= Ua.x*nx + Ua.y* ny + Ua.z* nz;
    double un = u_c * nx + v_c * ny + w_c * nz;
    double ut = u_c - un * nx, vt = v_c - un * ny, wt = w_c - un * nz;
    double ut_mag = sqrt( ut*ut + vt*vt + wt*wt );
    //*ustar = find_utau_Cabot_roughness(1./simCtx->ren, ut_mag, sc, 0.01, 0, ks);
        *ustar = find_utau_hydset(1./simCtx->ren, ut_mag, sc, 0.001, ks);
    //double ut_mag_modeled = u_Cabot_roughness(1./simCtx->ren, sb, *ustar, 0, ks);
    double ut_mag_modeled = u_hydset_roughness(1./simCtx->ren, sb, *ustar, ks);
    
    if (ut_mag_modeled<0.) {
      //      freeslip(user, sc,sb,Ua, Uc,Ub,nx,ny,nz);
      ut_mag_modeled=ut_mag;
      //*ustar =0.;
      //  LOG_ALLOW(LOCAL,LOG_DEBUG, "Y+>300!!!!!!!!!!!!!!! %le ",*ustar);
    }
    
    //PetscReal yp=sc* (*ustar)*simCtx->ren;
    //if (yp>300) LOG_ALLOW(LOCAL,LOG_DEBUG, "  Y+>300!!!!!!!!!!!!!!!  %le ",yp);
    
    if(ut_mag>1.e-10) {
        ut *= ut_mag_modeled/ut_mag;
        vt *= ut_mag_modeled/ut_mag;
        wt *= ut_mag_modeled/ut_mag;
    }
    else ut=vt=wt=0;
                    
    // u = ut + (u.n)n
    (*Ub).x = ut + (sb/sc * un ) * nx;
    (*Ub).y = vt + (sb/sc * un ) * ny;
    (*Ub).z = wt + (sb/sc * un ) * nz;
    
    (*Ub).x += Ua.x;
    (*Ub).y += Ua.y;
    (*Ub).z += Ua.z;
    
}

Here is the call graph for this function:

Here is the caller graph for this function:

wall_function_Cabot()

void wall_function_Cabot	(	UserCtx *	user,
		double	ks,
		double	sc,
		double	sb,
		Cmpnts	Ua,
		Cmpnts	Uc,
		Cmpnts *	Ub,
		PetscReal *	ustar,
		double	nx,
		double	ny,
		double	nz,
		double	dpdx,
		double	dpdy,
		double	dpdz,
		int	count
	)

Definition at line 319 of file wallfunction.c.

        {
            SimCtx *simCtx = user->simCtx;
        double u_c = Uc.x - Ua.x, v_c = Uc.y - Ua.y, w_c = Uc.z - Ua.z;
       //   double ua_n= Ua.x*nx + Ua.y* ny + Ua.z* nz;
            double un = u_c * nx + v_c * ny + w_c * nz;
            double ut = u_c - un * nx, vt = v_c - un * ny, wt = w_c - un * nz;
            double ut_mag = sqrt( ut*ut + vt*vt + wt*wt );         
            double dpdt = dpdx*ut/ut_mag + dpdy*vt/ut_mag + dpdz*wt/ut_mag;
            double dpdtn = sqrt(dpdx*dpdx+dpdy*dpdy+dpdz*dpdz-pow(dpdx*ut/ut_mag,2.0) - pow(dpdx*nx,2.0) - 
            pow(dpdy*vt/ut_mag,2.0) - pow(dpdy*ny,2.0) - pow(dpdz*nz,2.0)-pow(dpdz*wt/ut_mag,2.0));
            double utau,taw1,taw2;
            dpdtn=0.0; //dpdt=0.0;
               if (count ==0 ||count > 4){
            find_utau_Cabot(1./simCtx->ren, ut_mag, sc, 0.01, dpdt, dpdtn, &utau , &taw1, &taw2);
            *ustar=utau;}
            //*ustar=0.045;
            double ut_mag_modeled = u_Cabot(1./simCtx->ren, sb, *ustar, dpdt, taw1);
 
                    if(ut_mag>1.e-10) {
        ut *= ut_mag_modeled/ut_mag;
        vt *= ut_mag_modeled/ut_mag;
        wt *= ut_mag_modeled/ut_mag;
         }
          else ut=vt=wt=0;
                    
    // u = ut + (u.n)n
    (*Ub).x = ut + (sb/sc * un ) * nx;
    (*Ub).y = vt + (sb/sc * un ) * ny;
    (*Ub).z = wt + (sb/sc * un ) * nz;
    
    (*Ub).x += Ua.x;
    (*Ub).y += Ua.y;
    (*Ub).z += Ua.z;
                   
       }

Here is the call graph for this function: