CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > ANSYS > FLUENT > Fluent UDF and Scheme Programming

dynamic contact angle udf returns no value to solver

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   November 1, 2012, 11:32
Default dynamic contact angle udf returns no value to solver
  #1
New Member
 
Join Date: Aug 2011
Posts: 3
Rep Power: 15
shiraz_man67 is on a distinguished road
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)
}
shiraz_man67 is offline   Reply With Quote

Old   February 23, 2018, 10:45
Default
  #2
New Member
 
Join Date: Feb 2018
Posts: 6
Rep Power: 8
FelixJJ is on a distinguished road
Hi shiraz, have you solved your problem? I have a similar problem as you did.
FelixJJ is offline   Reply With Quote

Old   July 1, 2018, 04:33
Default
  #3
New Member
 
Join Date: Jul 2018
Posts: 7
Rep Power: 8
Albert Lee is on a distinguished road
Quote:
Originally Posted by shiraz_man67 View Post
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!
Albert Lee is offline   Reply With Quote

Old   July 1, 2018, 17:01
Default
  #4
New Member
 
Join Date: Feb 2018
Posts: 6
Rep Power: 8
FelixJJ is on a distinguished road
changing the unit of contact angle to rad in your code might help
FelixJJ is offline   Reply With Quote

Old   July 3, 2018, 03:11
Default
  #5
New Member
 
Join Date: Jul 2018
Posts: 7
Rep Power: 8
Albert Lee is on a distinguished road
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
Albert Lee is offline   Reply With Quote

Old   July 3, 2018, 14:51
Default
  #6
New Member
 
Join Date: Feb 2018
Posts: 6
Rep Power: 8
FelixJJ is on a distinguished road
yes, it's weird. And seems no one points it out clearly..... don't know why...
FelixJJ is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
dynamic contact angle implementation q doubtsincfd OpenFOAM Programming & Development 7 November 12, 2015 08:07
Slug Flow, interFoam, problems with Contact Angle PrzemekPL OpenFOAM Running, Solving & CFD 13 February 18, 2014 22:10
Dynamic Mesh solver OF_NACA OpenFOAM 0 June 21, 2011 03:03
Working directory via command line Luiz CFX 4 March 6, 2011 20:02
Dynamic contact angle with Fluent 6.3.26 baechtel FLUENT 0 May 21, 2009 23:55


All times are GMT -4. The time now is 20:43.