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

 shiraz_man67 November 1, 2012 12: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*/
{
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);
{
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);
}
{
real VOF_Val[3], VOF_Mag, source, VOF_Norm[3];
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. */
{
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;
}
}
source = 0.0;
}
{
real 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;
/* 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;
/* 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++;
x_Bisect=(x_Bottom+x_Top)/2.0;
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);
/* 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 12:16.