UDF for Dynamic contact angle
-
-
Hey I'm working on the simulation of the droplet impact on different surfaces
I need help with the udf code for the calculation of the dynamic contact angle
I'm using 2D axisymmetric model
I have written a udf code for that which is
#include "udf.h"
#include "mem.h"
/* Dynamic Contact UDF */
/* Enumerate UDM numbers */
enum{
VOF_G_X,
VOF_G_Y,
VOF_G_Z,
CONTACT_ANGLE,
CONTACT_N_REQ_UDM
};
#define VISCOSITY 0.001153
#define SURF_TENS 0.0727
#define static_Con_Ang 93.0
#define dyn_adv_Ang 108.0
#define dyn_rec_Ang 78.0
#define C_VOF_G_X(C,T)C_UDMI(C,T,VOF_G_X)
#define C_VOF_G_Y(C,T)C_UDMI(C,T,VOF_G_Y)
#define C_VOF_G_Z(C,T)C_UDMI(C,T,VOF_G_Z)
#define C_CONTACT_ANGLE(C,T)C_UDMI(C,T,CONTACT_ANGLE)
DEFINE_ADJUST(dynamic_contact, domain)
{
  Thread *ct, *pct, *t0;/*pointer to the thread*/
  Thread *ft;
  cell_t c, c0;/*An integer data type that identifies a particular cell within a cell thread.*/
  face_t f;/*An integer data type that identifies a particular face within a face thread*/
  double theta_n;
  real vel_re[ND_ND],Unit_tan[ND_ND], vof_n[ND_ND],VOF_Norm[ND_ND], vec_nw[ND_ND], vec_tw[ND_ND];
  real A[ND_ND], Amag, VOF_Mag,Tan_Mag, veltan[ND_ND];
  double f_Hoff_inverse, x_hoff, Ca;
  int phase_domain_index = 1;/*The phase_domain_index can be used in UDFs to distinguish between the primary and secondary
  phase-level threads. phase_domain_index is always assigned the value of 0 for the primary phaselevel thread*/
  Domain *pDomain = DOMAIN_SUB_DOMAIN(domain, phase_domain_index);
/*to get the specific phase (or subdomain) pointer within the mixture domain*/
  if (first_iteration)
  {
/*first iteration is true only at the first iteration of a timestep*/
/*Check UDMs have been set up */
    if (N_UDM < CONTACT_N_REQ_UDM)
    {
      Message0(" WARNING: Require at least %d UDMs to be set up for dynamic contact angle model. ", CONTACT_N_REQ_UDM);
      return;
    }
/*Calculate VOF gradient*/
    Alloc_Storage_Vars(pDomain,SV_VOF_RG,SV_VOF_G,SV_NULL); /* Primary storage of variables being calculated */
    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);
/* Store VOF gradient in user memory */
    thread_loop_c(ct, domain)
    {
 /* Simpler thread loop when you want to loop over all cell threads in a given domain */
      if (FLUID_THREAD_P(ct))
      {
        /*to check whether a cell thread is a fluid thread returns 1 (or TRUE) if the thread that is passed is a fluid thread, and 0 (or FALSE) if it is not*/
        pct = THREAD_SUB_THREAD(ct, phase_domain_index); /* Use this instead of pt [] */
        begin_c_loop(c, ct)
        {
          ND_V(C_VOF_G_X(c, ct), C_VOF_G_Y(c, ct),C_VOF_G_Z(c, ct ),=,C_VOF_G(c, pct)); /* Set 3 components to a vector */
        }
        end_c_loop(c, ct)
      }
/* Free memory used for VOF gradient */
      Free_Storage_Vars(pDomain,SV_VOF_RG,SV_VOF_G,SV_NULL);
    }
/* Calculate wall contact angle and save to user memory */
    thread_loop_f(ft, domain)
    {
/* to loop over all face threads in a given domain.*/
    if(THREAD_TYPE(ft) == THREAD_F_WALL)  /* Mixture Thread and Wall */
      {
        t0 = THREAD_T0(ft);
    /*used to identify cell threads that are adjacent to a given face f in a face thread t. THREAD_T0 expands to a function that returns the cell thread of a given face’s adjacent cell c0*/
        if (FLUID_THREAD_P(t0)) /* Fluid Zone */
          {
/* Set contact angle parameters according to the type of wall surface */
            begin_f_loop(f, ft)
            {
              c0 = F_C0(f, ft);
  /*used to identify cells that are adjacent to a given face thread t. F_C0 expands to a function that returns the index of a face’s neighboring c0 cell*/
/* Get flow velocity */
              NV_D(vel_re ,= ,C_U(c0, t0), C_V(c0, t0), C_W(c0, t0));
              F_AREA(A, f, ft) ; /* Area vector */
              Amag = NV_MAG(A);
              NV_S(A,/=,Amag);
              /*A now holds outward face normal*/
              NV_D(vof_n, =, C_VOF_G_X(c0, t0), C_VOF_G_Y(c0, t0 ), C_VOF_G_Z(c0, t0));
              VOF_Mag=NV_MAG(vof_n);
              if(VOF_Mag!=0.0)
              {
                VOF_Mag=NV_MAG(vof_n);
                NV_VS(VOF_Norm, =, vof_n, *, (1/VOF_Mag));
              }
              else
              {
                NV_VS(VOF_Norm, =, vof_n, *, (0.0));
              }
/* vof along wall*/
              NV_VS(vec_nw, =, A, *, NV_DOT( VOF_Norm, A)); /*normal relative velocity vector */
              NV_VV(vec_tw, =, VOF_Norm, -, vec_nw); /* Tangent velocity direction vector */
              Tan_Mag = NV_MAG(vec_tw);
              if(Tan_Mag!=0.0)
              {
                Tan_Mag = NV_MAG(vec_tw);
                NV_VS(Unit_tan, =, vec_tw, *, (1/Tan_Mag));
                NV_VS(veltan, =, Unit_tan, *, NV_DOT( vel_re, Unit_tan));
              }
              else
              {
                NV_D(veltan, =, C_VOF_G_X(c0, t0), C_VOF_G_Y(c0, t0 ), C_VOF_G_Z(c0, t0));
              }
/*Contact line direction*/
/*Receding : Calculate dynamic contact angle*/
              if (NV_DOT(veltan, vof_n) > 0.0)
                {
                  Ca = (VISCOSITY*(NV_MAG(veltan)))/SURF_TENS;
/*DRCA model*/
                  theta_n = ((pow((pow(((dyn_rec_Ang)*M_PI/180.0), 3) - 72.0*Ca),(1/3)))*180)/M_PI;
                  C_CONTACT_ANGLE(c0, t0) = theta_n ; /* deg */
                }
              else
                {
/* Advancing : Calculate dynamic contact angle*/
                  Ca = (VISCOSITY*(NV_MAG(veltan)))/SURF_TENS;
/* Kistler’s model*/
                  f_Hoff_inverse = pow(((dyn_adv_Ang)*M_PI/180.0), 3)/72.0;
                  x_hoff = Ca + f_Hoff_inverse;
                  theta_n = ((acos(1.0 - (2.0*tanh(5.16*pow((x_hoff/(1.0 + 1.31*pow(x_hoff, 0.99))),0.706)))))*180)/M_PI;
                  C_CONTACT_ANGLE(c0, t0) = theta_n; /* [ degs ] */
                }
            }
            end_f_loop (f, ft)
          }
      }
    }
  }
  printf("Executed Adjust ");
}
DEFINE_PROFILE(contact_angle, t, i)
{
  face_t f;
  cell_t c0;
  Thread *t0;
  begin_f_loop(f, t)
    {
    c0 = F_C0(f ,t);
    t0 = F_C0_THREAD(f , t );
    F_PROFILE(f ,t ,i) = C_CONTACT_ANGLE(c0, t0); /* Angle within the liquid in degrees */
    }
  end_f_loop(f, t)
}
DEFINE_EXECUTE_AFTER_CASE(set_name, libname )
{
Message0 (" Setting UDM names . . . ");
  Set_User_Memory_Name(VOF_G_X,"VOF Gradient−x");
  Set_User_Memory_Name(VOF_G_Y,"VOF Gradient−y");
  Set_User_Memory_Name(VOF_G_Z,"VOF Gradient−z ");
Set_User_Memory_Name(CONTACT_ANGLE," Contact Angle [deg]");
  Message0 ("Done. ") ;
}
But I'm getting an error after loading the udf and assigning the DEFINE_ PROFILE udf
the error is following
Node 0: Process 23624: Received signal SIGSEGV.
================================================== ============================
999999: mpt_accept: error: accept failed: No error
999999: mpt_accept: error: accept failed: No error
999999: mpt_accept: error: accept failed: No error
999999: mpt_accept: error: accept failed: No error
999999: mpt_accept: error: accept failed: No error
999999: mpt_accept: error: accept failed: No error
.
.
.
.
.
.
.
The fl process could not be started.
Please help me out, I'm new with these udfs
Viewing 4 reply threads- The topic ‘UDF for Dynamic contact angle’ is closed to new replies.
-
You are navigating away from the AIS Discovery experience
