CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > 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 Display Modes
Old   November 1, 2012, 12:32
Default dynamic contact angle udf returns no value to solver
  #1
New Member
 
Join Date: Aug 2011
Posts: 3
Rep Power: 5
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

Reply

Thread Tools
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 On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
dynamic contact angle implementation q doubtsincfd OpenFOAM 6 May 5, 2015 09:30
Slug Flow, interFoam, problems with Contact Angle PrzemekPL OpenFOAM Running, Solving & CFD 13 February 18, 2014 23:10
Dynamic Mesh solver OF_NACA OpenFOAM 0 June 21, 2011 03:03
Working directory via command line Luiz CFX 4 March 6, 2011 21: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 15:07.