[shiraz_man67]

I wrote dynamic contact angle for wall in VOF method but the function most of times returns no value to the solver . 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)
}

[FelixJJ]

Hi shiraz, have you solved your problem? I have a similar problem as you did.

[Albert Lee]

Originally Posted by shiraz_man67

I wrote dynamic contact angle for wall in VOF method but the function most of times returns no value to the solver . 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)
}

Hi shiraz, have you solved your problem? What I studied was the capillary phenomenon in a tube partially filled with water. When I builded and loaded a udf of dynamic contact angle, the setting process showed no errors. But the contact angle in iteration was not what I setted, even though the simplest udf file as you posted at the last of your question.

I have being stuck in this problem for several days! Nothing help could come out!

If you have got the correct method for setting for coding the udf of dynamic contact angle, please give me a reply.

Thank you very much!

[FelixJJ]

changing the unit of contact angle to rad in your code might help

[Albert Lee]

It really works! But I don't know why this could happen. And FLUENT didn't illustrate it in help files. Maybe it's a bug？ Anyway，
Thank you very much!

Albert Lee

[FelixJJ]

yes, it's weird. And seems no one points it out clearly..... don't know why...