关于VOF中调用volume fraction梯度的问题

nictse911

大家好，我想用UDF把一个新的蒸发相变模型加入VOF模型中来计算液滴蒸发，液滴直径为100微米，初始温度为500K，周围空气初始温度为800K。网格尺寸1微米，步长0.1微秒。

公式中需要用到volume fraction的梯度。我用的C_VOF_G(cell, gas)，编译能通过，也可以进行计算。但是计算结果中mass transfer rate 始终为0，然后我加入了printf来查看公式中的中间值，发现用C_VOF_G(cell, gas)获取的界面上的梯度值始终为0，而后处理中界面上的VOF gradient不是0。

UDF如下：

1. /* UDF to define a simple mass transfer based on SunDongliang's equation. */
   
2.    
   
3.   #include "udf.h"
   
4.   #include "sg_mphase.h"
   
5. 
   
6.   DEFINE_MASS_TRANSFER(liq_gas_source_Sun, cell, thread, from_index, from_species_index, to_index, to_species_index)
   
7.   {
   
8.     real m_lg; // mass transfer rate on the interface, kg/s/m3 ;
   
9.     real T_SAT = 373.15; // saturated temperature, K
   
10.     real h_lg = 2257600; // evaporation enthalpy, J/kg
    
11.     real vof_cutoff = 0.05;
    
12.     real vpVofT; // dot product of temperature gradient and volume fraction gradient
    
13.     real v_sx; // volume fraction gradient dx
    
14.     real v_sy; // volume fraction gradient dy
    
15.     Thread *liq, *gas;   
    
16.     liq = THREAD_SUB_THREAD(thread, from_index);
    
17.     gas = THREAD_SUB_THREAD(thread, to_index); //  from_index and to_index are defined in the Mass panel of Phase interaction.
    
18.     m_lg = 0.0; // Initial value of m_lg.
    
19.    
    
20. 
    
21.     if (NULL != THREAD_STORAGE(gas,SV_VOF_G))
    
22.     {
    
23.      v_sx = C_VOF_G(cell, gas)[0];
    
24.      v_sy = C_VOF_G(cell, gas)[1];     
    
25.      }
    
26.      
    
27.     if ((C_VOF(cell, liq) > vof_cutoff) && (C_VOF(cell, liq) < (1-vof_cutoff)))
    
28.         {                                               // Evaporating, and it will be added to vapor phase governing equation
    
29.         vpVofT = v_sx*C_T_G(cell, thread)[0] + v_sy*C_T_G(cell, thread)[1];
    
30.           printf("v_sy: %f \n", v_sy);
    
31.           printf("TempGradx: %f \n", C_T_G(cell, thread)[0]);
    
32.           printf("dot product of TempGrad and VOFGrad: %f \n", vpVofT);
    
33.         m_lg = 2/h_lg*C_K_L(cell, gas)*vpVofT; //NV_DOT(C_VOF_G(cell, gas), C_T_G(cell, thread));
    
34.         }
    
35.        return (m_lg);
    
36.    }
    
37. 
    

复制代码

麻烦各位给分析下，先行谢过！