[HFLUENT]

Hi,

I try to simulate wall slipping of a polymer melt by using Fluent. The idea is that the melt starts slipping at the wall above a critical shear stress (Tau_krit). Above this cirtical shear stress the wall shear stress follows the equation

Tau = Tau_krit + K*U_gleit,

where K is an material parameter (slipping coefficient) and U_gleit is the slipping velocity. For this im using the UDF attached below.

My problem is, that the equation mentioned above is only valid for the case if Tau>Tau_krit, otherwise the "no slip" condition of Fluent should be used for the calculation. Is there a possibility to affect this setting by an UDF? Up to know I am using the following expression to regard the mentioned condition in my UDF:

if (taumag111 < tau_krit0) //taumag111=calculated shear stress
F_PROFILE(f,t,i)= tau0[dir]; //tau0=shear stress in case of wall adhesion
else
F_PROFILE(f,t,i) = tau[dir]; //tau=shear stress in case of wall slipping

But this does not provide a useful result of my simulation. I think the Problem is the specified shear boundary condition affected by this approach.
Another Idea is to set the velocity of the cells close to the boeundary to zero. Would that work? How can I include this in the UDF?

It would be great if someone would have an Idea.

Thank You!


To get an overview, the whole UDF:

/* Calculation of wall slipping following: */
/* Tau = Tau_krit + K*U_gleit */

/* Fluent-Header */
#include "udf.h"
#include "surf.h"

void shear_profile(Thread* t, int i, int dir)
{
face_t f; //face t is the data type for a face ID. It is an integer index that identifes a particular face within a given thread.
Thread* t0 = THREAD_T0(t); //The Thread data type is a structure that acts as a container for data that is common to the group of cells or faces that it represents. For multiphase applications, there is a thread structure for each phase, as well as for the mixture. See Section 1.8.1 for details.
real omega;
real mu_c;
real S_dot;
real w_z;
real NV_VEC(du); /* declaring vector du */
/*******************************************/
#if !RP_NODE
omega = RP_Get_Real("mom/relax");
#endif
host_to_node_real_1(omega);
/*******************************************/

/* loop on all faces on the wall */
begin_f_loop (f,t)
{
cell_t c0 = F_C0(f,t); /* neighboring cells of face f, and their (corresponding) threads */
real A[ND_ND], rc[ND_ND], rf[ND_ND], dr[ND_ND], vc[ND_ND],vcf[ND_ND], At, delta; // The constant ND ND is defined as 2 for RP 2D (FLUENT 2D) and RP 3D (FLUENT 3D). It can be used when you want to build a 2 x 2 matrix in 2D and a 3 x 3 matrix in 3D. When you use ND ND, your UDF will work for both 2D and 3D cases, without requiring any modifications.
real NV_VEC(deltavec); /* declaring vector deltavec */
real NV_VEC(vc_wall); /* declaring vector vc_wall */
real NV_VEC(wall_normal_unit_vect); /* declaring vector wall_normal_unit_vect */
real NV_VEC(vc_wall_normaltowall); /* declaring vector vc_wall_normaltowall */
real NV_VEC(vc_wall_projected); /* declaring vector vc_wall_projected */
real mytmp;
real NV_VEC(vc_wall_projected_UNIT); /* declaring vector vc_wall_projected_UNIT */
real tau_krit0 = 1000000000000000000;
real k_konstant = 10;
real NV_VEC(vc_wall_projected_voll); /* declaring vector vc_wall_projected_voll */
real NV_VEC(vc_wall_projected_tau); /* declaring vector vc_wall_projected_tau */
real NV_VEC(S_dot1); /* declaring vector S_dot1 */
real NV_VEC(tau); /* declaring vector tau */
real NV_VEC(tau0); /* declaring vector tau0 */
real taumag111=NV_MAG(tau0); //The utility NV MAG2 computes the sum of squares of vector components.
real NV_VEC(temp_for_wall);

/* Get normal distance between face and neighboring cell. */
F_AREA(A,f,t); /* returns areavector; Ax, Ay, Az */
At = NV_MAG(A); /* Area magnitude / Flaechenbetrag */
F_CENTROID(rf,f,t); /* face Centroid; xf, yf, zf */
C_CENTROID(rc,c0,t0); /* cell centroid; xc, yc, zc*/
NV_VV(dr,=,rf,-,rc); /* distance between face and cell centroid, (vector) The utility NV_VV performs operations on vector elements. The operation that is performed on the elements depends upon what symbol ( -,/,*) is used as an argument in place of the + signs in the following macro call. */
delta = fabs(NV_DOT(dr,A)/At); /* normal distance (dr*A*cos phi ) */

/*================================================= =================*/
// V denotes a vector, S denotes a scalar, and D denotes a sequence of three vector components of which the third is always ignored for a two-dimensional calculation
NV_V(deltavec, =, A); /* set deltavec equal A (vectorial) The utility NV_V performs an operation on two vectors. */
NV_S(deltavec, /=, At); /* get unit area vector*/
NV_S(deltavec, *=, delta); /* get vector from cell center to wall*/

/*Calculation of the velocity at the wall using velocity vector vc */

NV_D(vc,=,C_U(c0,t0),C_V(c0,t0),C_W(c0,t0)); /* defining vc in terms of velocity field ?? */

/*Defining an wall velocity
NV_D(vc, +=, 0,0,0/*temp_for_wall*/);
/* ************************************************** ********** */

NV_V(vc_wall, =, vc); /* vc_wall = vc */
NV_D(vc_wall, +=, NV_DOT(C_U_G(c0, t0), deltavec), /* dot product of the two returned */
NV_DOT(C_V_G(c0, t0), deltavec),
NV_DOT(C_W_G(c0, t0), deltavec));
NV_D(vcf,=,F_U(f,t),F_V(f,t),F_W(f,t));

/* Calculation of the dirction vector of critical wall shear stress \tau_kr
* The velocity is already in the vairable vc, the velocity extrapolated to the wall vc_wall.
* This has to projected to the wall surface. Hence the component normal to the wall is determined and subtract this */
NV_VS(wall_normal_unit_vect, =, A, /, At); // /* Normalisieren?? */
mytmp = NV_DOT(vc_wall, wall_normal_unit_vect); /* hier ! */
NV_VS(vc_wall_normaltowall, =, wall_normal_unit_vect, * , mytmp);
NV_VV(vc_wall_projected, =, vc_wall, -, vc_wall_normaltowall);
NV_VV(du, =, vc_wall_projected, -, vcf);

/* Now a direction vector of length=1 is generated from the projected velocity*/
mytmp = NV_MAG(du);
NV_VS(vc_wall_projected_UNIT, =, du, /, mytmp);

/* This is the direction vector (length: 1) for the cirtical wall shear stress.
* Multiplication of the absolute value of wall shear stress
* with the x,y,z-direction-components of "vc_wall_projected_UNIT" to get the x,y,z-direction-components
* of the critical wall shear stress \tau_kr */

/*here the implemented equation Tau = Tau_krit + K. U_gleit */
NV_VS(vc_wall_projected_tau, =, vc_wall_projected_UNIT, *, tau_krit0); // vc_wall_prj_tau = vc_wall_projected_UNIT * tau_krit0 /*Tau_krit0 projected to the velocity components */

NV_VS(vc_wall_projected_voll, =, du, *, k_konstant); /**/

NV_VV(tau, =, vc_wall_projected_tau, +, vc_wall_projected_voll) ; /**/

/* if slipping dosent occur, calculation with FLUENT shear stress values */
if (NNULLP(THREAD_STORAGE(t, SV_WALL_SHEAR))) // If SV_WALL_SHEAR is allocated
{
/* Calculation of the laminar shear stress in the three directions S_dot*/
mu_c = C_MU_L(c0,t0); /* retrieve cell laminar viscosity */
S_dot = C_STRAIN_RATE_MAG(c0,t0); //strain/shear rate magnitude
NV_VS(S_dot1, =, vc_wall_projected_UNIT, *, S_dot); // ???? S_dot1 wird nicht zugewiesen und verwendet // Skalarmultiplikation

/*Calculation of shear stress by velocity components and wall distance*/
tau0[0] = mu_c * du[0]/delta ;
tau0[1] = mu_c * du[1]/delta;
tau0[2] = mu_c * du[2]/delta;
}
else
{
printf("Error\n");
}

/*F PROFILE can be used to store a boundary condition in memory for a given face and
thread, and is typically nested within a face loop. See mem.h for the complete macro
definition for F PROFILE.
*/
if (taumag111 < tau_krit0)
F_PROFILE(f,t,i)= tau0[dir];
else
F_PROFILE(f,t,i) = tau[dir];
}
end_f_loop(f,t)
}

// DEFINE_MACRONAME(udf_name, passed-in variables)
// e.g: DEFINE_PROFILE(inlet_x_velocity, thread, index)
DEFINE_PROFILE(x_shear_profile, t, i)
{
shear_profile(t, i, 0);
}


DEFINE_PROFILE(y_shear_profile, t, i)
{
shear_profile(t, i, 1);
}


DEFINE_PROFILE(z_shear_profile, t, i)
{
shear_profile(t, i, 2);
}