CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   Fluent UDF and Scheme Programming (http://www.cfd-online.com/Forums/fluent-udf/)
-   -   dynamic contact angle udf returns no value to solver (http://www.cfd-online.com/Forums/fluent-udf/108792-dynamic-contact-angle-udf-returns-no-value-solver.html)

shiraz_man67 November 1, 2012 13:32

dynamic contact angle udf returns no value to solver
 
I wrote dynamic contact angle for wall in VOF method but the function most of times returns no value to the solver :confused: . i even try to write simplest case of udf az follow bur the contact angle I see is still wrong here the complete code:
#include "udf.h"
#include "sg_mphase.h"
#include "sg_vof.h"
#include "sg.h"
#include "mem.h"
#define VISCOSITY 0.001
#define SURF_TENS 0.0728
#define MYTRUE 1
#define MYFALSE 0
#define Hoff(x) acos(1-(2.0*tanh(5.16*(pow((x/(1+(1.31*pow(x, 0.99)))), 0.706)))))
#define static_Con_Ang 30.
#define index_source 3
/* This Code computes the normals of the VOF function*/
DEFINE_ADJUST(store_gradient, domain)
{
Thread *t;
Thread **pt;
cell_t c;
int phase_domain_index = 1; /* Secondary Domain */
Domain *pDomain = DOMAIN_SUB_DOMAIN(domain,phase_domain_index);
void calc_source();
Alloc_Storage_Vars(pDomain,SV_VOF_RG,SV_VOF_G,SV_N ULL);
Scalar_Reconstruction(pDomain, SV_VOF,-1,SV_VOF_RG,NULL);
Scalar_Derivatives(pDomain,SV_VOF,-1,SV_VOF_G,SV_VOF_RG,
Vof_Deriv_Accumulate);
mp_thread_loop_c (t,domain,pt)
if (FLUID_THREAD_P(t))
{
Thread *ppt = pt[phase_domain_index];
begin_c_loop (c,t)
{
calc_source(c,t);
}
end_c_loop (c,t)
}
Free_Storage_Vars(pDomain,SV_VOF_RG,SV_VOF_G,SV_NU LL);
}
void calc_source(cell_t cell, Thread *thread)
{
real VOF_Val[3], VOF_Mag, source, VOF_Norm[3];
Thread *phaset;
phaset= THREAD_SUB_THREAD(thread,1);
if(C_VOF(cell,phaset)!=0.0 && N_TIME > 1)
{
/* The gradients of the VOF function are found in the x,y and z dir. */
if (NULLP(THREAD_STORAGE(phaset, SV_VOF_G)))
{
Message0("N_TIME = %d, ....show-grad:Gradient of VOF is not available \n ", N_TIME);
Error("0");
}
VOF_Val[0]=-C_VOF_G(cell,phaset)[0];
VOF_Val[1]=-C_VOF_G(cell,phaset)[1];
VOF_Val[2]=-C_VOF_G(cell,phaset)[2];
/* The magnitude of the VOF gradients is found so it can be normalized */
VOF_Mag=NV_MAG(VOF_Val);
if(VOF_Mag!=0.0)
{
VOF_Mag=NV_MAG(VOF_Val);
VOF_Norm[0]=VOF_Val[0]/VOF_Mag;
VOF_Norm[1]=VOF_Val[1]/VOF_Mag;
VOF_Norm[2]=0;
}
else
{
/* This is to avoid the divide by zero function*/
VOF_Norm[0]=0.0;
VOF_Norm[1]=0.0;
VOF_Norm[2]=0.0;
}
C_UDMI(cell,thread,0)=VOF_Norm[0];
C_UDMI(cell,thread,1)=VOF_Norm[1];
C_UDMI(cell,thread,2)=0;
}
source = 0.0;
C_UDMI(cell, thread, index_source) = source;
}
DEFINE_SOURCE(VOF_Norms, cell, thread, dS, eqn)
{
real source;
source = C_UDMI(cell, thread, index_source);
return source;
}
/* This Define_profile code is designed to provide a dynamic contact angle for the VOF function*/
DEFINE_PROFILE(con_ang, t, pos)
{
/* First the various pointer variables are created*/
face_t f;
cell_t c;
real feta_d, vel_Val[3], cont_Line_Vel, VOF_Normal[3], cap_Num, static_Con_Rad, x_Bottom,h,
x_Top, x_Bisect, hoff_Old, hoff_Cur, hoff_New, finish_Cond, inv_Hoff=0.0;
int notConverged, itNum;
Thread *t0,*pt;
/* This code is designed to find the zero for the inverted hoffman function by finding the zero */
/* of the function at which the hoffman function results in the static contact angle */
/* First the variables are assigned*/
notConverged=MYTRUE;
x_Bottom=0.001;
x_Top=2.0;
itNum=0;
static_Con_Rad=((static_Con_Ang*M_PI)/180.);
/* A while loop performs the bisection method, a simple but very stable zero finder */
while(notConverged)
{
/* The variables used in the bisection method are assigned and the hoffman */
/* functions are evaluated */
itNum++;
hoff_Old=(Hoff(x_Bottom)- static_Con_Rad);
hoff_Cur=(Hoff(x_Top)- static_Con_Rad);
x_Bisect=(x_Bottom+x_Top)/2.0;
hoff_New=(Hoff(x_Bisect)- static_Con_Rad);
finish_Cond=fabs(1-(x_Bisect/x_Top));
/* The loop ends when the relative error is less than 1e-8 and the inverse */
/* hoffman value is stored for use later */
if(finish_Cond<0.00000001 || itNum>10000000)
{
inv_Hoff=x_Bisect;
notConverged=MYFALSE;
}
/* Conditions for the bisection method */
if((hoff_Old*hoff_New)<0.0)
{
x_Top=x_Bisect;
}
if((hoff_Cur*hoff_New)<=0.0)
{
x_Bottom=x_Bisect;
}
}
/* Now the main loop goes through all the faces in the boundary */
begin_f_loop(f,t)
{
/* The cell and phase threads are isolated */
c=F_C0(f,t);
t0=THREAD_T0(t);
pt= THREAD_SUB_THREAD(t0,1);
/* The main formulation is only applied if the VOF is >0 */
if(C_VOF(c,pt)!=0.0)
{
/* The velocities are recorded in each direction */
vel_Val[0]=C_U(c,t0);
vel_Val[1]=C_V(c,t0);
vel_Val[2]=0;
/* The VOF normals are brought in */
VOF_Normal[0]=C_UDMI(c,t0,0);
VOF_Normal[1]=C_UDMI(c,t0,1);
VOF_Normal[2]=0;
/* The contact line vel. is calc from the dot product of VOF and Vel */
cont_Line_Vel=NV_DOT(vel_Val,VOF_Normal);
/* The capillary number is found based on cont line vel. */
cap_Num=((VISCOSITY*cont_Line_Vel)/SURF_TENS);
/* The dynamic contact angle is defined then stored in the profile */
C_UDMI(c,t0,4)=cont_Line_Vel;

if ( cont_Line_Vel>0.000001) {
h=48;
}
else if (cont_Line_Vel<-0.000001 ) {
h=10;
}
else {
h=30;
}
C_UDMI(c,t0,5)=h;

feta_d=h;
F_PROFILE(f,t,pos)=feta_d;
}
else
{
F_PROFILE(f,t,pos)=static_Con_Ang;
}
C_UDMI(c,t0,6)=F_PROFILE(f,t,pos);
}
end_f_loop(f,t)
}

and the simplest case code with just a constant value:
#include "udf.h"
DEFINE_PROFILE(CA_profile,t,i)
{
face_t f;
begin_f_loop_all(f,t)
{
F_PROFILE(f,t,i)=30.0;
}
end_f_loop_all(f,t)
}


All times are GMT -4. The time now is 22:53.