CFD Online Discussion Forums - wall slip law for newtonian fluid

Hi,

I am trying to input a wall slip law for a newtonian fluid in a pipeline.
The law has a yield stress limit $\tau_0$ .
When the yield stress is less than $\tau_0$ , there is adhesion of the fluid to the wall.
When the yield stress is greater than $\tau_0$ , there is a slip velocity at the wall like :
$u_{wall}=\mu*\frac{\partial u_z}{\partial r}$

The udf is here :

Code:

#define Cf 1                        /* friction coefficient */

#define TAU0 0.2                        /* limit yield stress for slipping */



DEFINE_PROFILE(wall_velocity_x, thread, position)

{

        face_t f;

        Thread *t0;

        cell_t c0;

        double VC[ND_ND];

        double TAUxy, VC_MAG, vg, mu;

        double dudy;

        begin_f_loop(f, thread)

                {

                t0 = THREAD_T0(thread);                /* adjacent cell thread to f */

                c0 = F_C0(f, thread); 

                dudy = C_DUDY(c0,t0);

                mu = C_MU_L(c0,t0);        

                VC[0]=-mu*dudy;

                VC[1]=0;                                /* stress vector's y componant*/

                if (NV_MAG(VC)==0)                /* stress vector's magnitude */

                        {VC_MAG = 0.00000001;}

                else {VC_MAG=NV_MAG(VC);}

                vg=VC[0]/Cf;

                if ( VC_MAG < TAU0)

                        {

                        F_PROFILE(f, thread, position)=0;

                        }

                else        {

                        F_PROFILE(f, thread, position)=vg;

                        }

                }

        end_f_loop(f, thread)

}

At first I run Fluent with a "wall" boundary condition. And after some steps, I change the wall to a "velocity inlet" boundary condition in order to put the udf as a component of the axial velocity.

But, the continuity residual goes to $10^{20}$ .

I think the reason is the fact that I use the derivative of the adjacent cell of the boundary to put it as the value of the face of the boundary. But I don't know how to make it in an other way.

Help me please !