动态接触角问题-有做VOF模拟自由液面的朋友一起来讨论下

liujun10607

我现在在做一个封闭环境中的流体流动状态的仿真，气液两项流，在仿真的时候定义一个region，然后江这个区域patch成液体，通过仿真来讨论液体在重力和表面张力的作用下形状随时间的变化，通过阅读文献我发现fluent6.3里面设置是接触角是静态接触角，而现在的一些研究多采用导入动态接触角的方式来完善fluent模型，下面是一篇最新的论文中的一接触角的程序，3D的我粘在下面和大家分享下：
#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 60.
#define index_source 2


/* This Code computes the normals of the VOF function*/

DEFINE_ADJUST(store_gradient, domain)
{
      Thread *t;
      Thread **pt;
          Thread *phaset;

      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_NULL);
      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_NULL);
      

}

void calc_source(cell_t cell, Thread *thread)
   {
   
  real VOF_Val[2], VOF_Mag, source, VOF_Norm[2];
  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];
  

    /* 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;
   
       }
    else
       {
      /* This is to avoid the divide by zero function*/
      VOF_Norm[0]=0.0;
      VOF_Norm[1]=0.0;
   
       }
    C_UDMI(cell,thread,0)=VOF_Norm[0];
    C_UDMI(cell,thread,1)=VOF_Norm[1];
   }
     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;

}

DEFINE_PROFILE(con_ang, t, pos)
{
  /* First the various pointer variables are created*/
  face_t f;
  cell_t c;
  real feta_d, vel_Val[2], cont_Line_Vel, VOF_Normal[2], cap_Num, static_Con_Rad, x_Bottom,
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);
      
      /* The VOF normals are brought in         
  */
      VOF_Normal[0]=C_UDMI(c,t0,0);
      VOF_Normal[1]=C_UDMI(c,t0,1);
      
      /* 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=fabs((VISCOSITY*cont_Line_Vel)/SURF_TENS);

      /* The dynamic contact angle is defined then stored in the profile  */
               if(cap_Num+inv_Hoff<0.0)
         {
        cap_Num=inv_Hoff;
         }
      feta_d=((Hoff(cap_Num+inv_Hoff))*180)/M_PI;
      F_PROFILE(f,t,pos)=feta_d;
  
      }
    else
      {
      F_PROFILE(f,t,pos)=static_Con_Ang;
      }
         }
  end_f_loop(f,t)
}

这个程序是能够编译成功的，也能够成功的导入到fluent6.3中，可是在迭代的时候报错，不知道是不是我udf下载的地方不对，还是什么？希望有做这方面研究的朋友能帮助我下，也希望这个程序对你们能有帮助~~

grmouse

请问这是哪篇文章中的代码？真心求教啊~

liujun10607

Liquid Water Dynamics in a Model Polymer Electrolyte Fuel Cell Flow Channel
上面讲的还蛮清晰的，程序是可以导进去可是好像那个接触角和Fluent中本身自带的表面张力模型（csf模型）有冲突，你仿成功了，可以教教我，呵呵~~

ndzhaolin

我也是做这方面的，必须要用动态接触角吗？

wesily

我现在在做水下气泡的运动模拟，有谁知道怎样生成随机气泡吗？

jiabinguang

霍夫曼的模型似乎还有接触角的限制，低于30度以下不能做