更圆润的液滴——UDF初始化

1 问题描述

通常，使用UDF进行VOF初始化时，假定网格中心位于球面（液滴）以内，则该网格充满液体，否则充满空气。该处理对于远离球面的网格是可行的，对于与球面相交的网格是不合适的。

图1给出了该处理的结果。从图可以看到，球面边缘非常粗糙。当网格进一步细化时，该问题可能会有缓解，但由于算法缺陷，该问题不可能根除。

图 1 一般UDF初始化结果

针对该问题，本文编制了一套算法，准确计算与球面相交网格内的液相体积分数，从而根除上述问题，结果如图2所示。

图 2 本文UDF初始化结果

2 一般初始化UDF

一般初始化UDF源代码：

#include "udf.h"

#include"math.h"

double Curve_Function(double *x)

{

/* 输入曲面方程，根据返回值判断点是否在曲面内部*/

double temp;

double a=3.0;  /* x 轴 半径*/

double b=3.0;  /* y 轴 半径*/

double c=3.0;  /* z 轴 半径*/

double d=1.0;  /* 方程常数 */

temp=pow((x[0]-0.0)/a,2.0)+pow((x[1]-0.0)/b,2.0)+pow((x[2]-0.0)/c,2.0)-d;

return temp;

} 

DEFINE_INIT(My_init,d)

{

   cell_t c;

   Thread *t;

   double X[3]; 

   double temp;

   thread_loop_c(t,d)

   {

        begin_c_loop(c,t)

        {

     C_CENTROID(X,c,t);

temp=Curve_Function(X);

if(temp>0.0)

{

C_UDMI(c,t,0)=0.0;

}

else

{

C_UDMI(c,t,0)=1.0; 

}

        }

        end_c_loop(c,t)

   }

} 

3 本文UDF

本文研制UDF源代码

#include "udf.h"

#include"stdlib.h"

#include"stdio.h"

#include"math.h"

#include"time.h"

double Curve_Function(double *x)

{

/* 输入曲面方程，根据返回值判断点是否在曲面内部*/

double temp;

double a=3.0;  /* x 轴 半径*/

double b=3.0;  /* y 轴 半径*/

double c=3.0;  /* z 轴 半径*/

double d=1.0;  /* 方程常数 */

temp=pow((x[0]-0.0)/a,2.0)+pow((x[1]-0.0)/b,2.0)+pow((x[2]-0.0)/c,2.0)-d;

return temp;

} 

double Get_Node_Coordinate(int c,Thread *t)

{

    Node *node;

int n;

int Flag1=0;

int Flag2=0;

int i;

double temp=0.0;

    double alpha;

double Point_Total=0;

double Point_In=0;

double Error=0.00001;

double Dev=1.0;

double Y[3]={0.0,0.0,0.0};

    double X[6]; /*0、1、2依次记录三个方向的最小值，3、4、5依次记录三个方向的最大值 */

c_node_loop(c,t,n)

{

node = C_NODE(c,t,n);

Y[0]=NODE_X(node);

Y[1]=NODE_Y(node);

Y[2]=NODE_Z(node);

temp=Curve_Function(Y);

if(temp<0.0)

Flag1=-1;

else if(temp>0.0)

Flag2=1;

}

/* 1, 判断顶点都在曲面内、曲面外或位于两侧*/

if((Flag1<-0.5)&&(Flag2<0.5))      /* 顶点全部在曲面内侧*/

{

alpha=1.0;

}

else if((Flag1>-0.5)&&(Flag2>0.5)) /*顶点全部在曲面外侧*/

{

alpha=0.0;

}

else if((Flag1<-0.5)&&(Flag2>0.5))/* 网格顶点位于 曲面两侧 */

{

/* 获取网格六个顶点坐标 */  

/* 2, 对于位于曲面两侧情形，获取网格坐标的范围，进而随机抽取点P*/

/* 2.1 获取网格顶点坐标,获取X、Y、Z三个方向的最大值最小值*/

        c_node_loop(c,t,n)

    {

    node = C_NODE(c,t,n);

    if(n==0)

    {

    X[0]=NODE_X(node);  X[0+3]=NODE_X(node);

X[1]=NODE_Y(node);  X[1+3]=NODE_Y(node);

X[2]=NODE_Z(node);  X[2+3]=NODE_Z(node);   

    }

    else

    {

    if(X[0]>NODE_X(node))X[0]=NODE_X(node); else if(X[0+3]<NODE_X(node)) X[0+3]=NODE_X(node);

    if(X[1]>NODE_Y(node))X[1]=NODE_Y(node); else if(X[1+3]<NODE_Y(node)) X[1+3]=NODE_Y(node);                               

    if(X[2]>NODE_Z(node))X[2]=NODE_Z(node); else if(X[2+3]<NODE_Z(node)) X[2+3]=NODE_Z(node);                          

    }

    }

while(Dev>Error)

{

    /* 随机抽取一个点 */

for(i=0;i<100;i++)

{

       // srand((unsigned)time(NULL));   /* Windows 下不适用 */

           temp=rand()/(RAND_MAX + 1.0); /* 0~1 之间的数     */

   // Message("temp =%lf\n",temp);

       Y[0]=X[0]*(1.0-temp)+X[0+3]*temp;

              temp=rand()/(RAND_MAX + 1.0); /* 0~1 之间的数 */

       Y[1]=X[1]*(1.0-temp)+X[1+3]*temp;

           temp=rand()/(RAND_MAX + 1.0); /* 0~1 之间的数 */

       Y[2]=X[2]*(1.0-temp)+X[2+3]*temp;

   temp=Curve_Function(Y);

               if(temp<0.0)

                 Point_In+=1;

               Point_Total+=1;   

}

/* 计算平均值  */

alpha= Point_In/Point_Total;

/* 计算标准偏差 */

temp=Point_In*(1.0-alpha)*(1.0-alpha)+(Point_Total-Point_In)*(alpha)*(alpha);

temp=temp/Point_Total;

Dev=sqrt(temp);

/* 强制退出 */

            if(Point_Total>1000000.0) break;

}

}

return alpha;

}

DEFINE_INIT(My_init,d)

{

   cell_t c;

   Thread *t;

   thread_loop_c(t,d)

   {

        begin_c_loop(c,t)

        {

C_UDMI(c,t,0)=Get_Node_Coordinate(c,t);

        }

        end_c_loop(c,t)

   }

} 

本文完。