fluent 曳力系数UDF



  • 液固两相双流体模型中,我写了一个Huilin-Giaspow曳力系数UDF,如下:

    #include "udf.h"
    #include "mem.h"
    #include "sg_mphase.h"
    #include "stdio.h"
    
    #define PAI 3.14159
    
    DEFINE_EXCHANGE_PROPERTY(Huilin_Giaspow,cell,mix_thread,s_col,f_col)
    {
    	real k_wenyu, k_ergun, k_l_s;
    	real x_vel_l, x_vel_s, y_vel_l, y_vel_s, z_vel_l, z_vel_s, vel_l, vel_s, abs_v;
    	real den_l, vis_l, vol_s;
    	real rey, diam_s, dcef, stf;
    	real g0, g1, g2, g3, g4;/*intermediate variable*/
    	
    	Thread *thread_l, *thread_s;
    	
    	/*find the threads for the liquid (primary)and solids(secondary phases)*/
    	thread_l = THREAD_SUB_THREAD(mix_thread, s_col);/*liquid phase*/
    	thread_s = THREAD_SUB_THREAD(mix_thread, f_col);/*solid phase*/
    	
    	/*find phase velocities*/
    	x_vel_l = C_U(cell, thread_l);
    	y_vel_l = C_V(cell, thread_l);
    	z_vel_l = C_W(cell, thread_l);
    	x_vel_s = C_U(cell, thread_s);
    	y_vel_s = C_V(cell, thread_s);
    	z_vel_s = C_W(cell, thread_s);
    	vel_l = sqrt(pow(x_vel_l,2.0)+pow(y_vel_l,2.0)+pow(z_vel_l,2.0));
    	vel_s = sqrt(pow(x_vel_s,2.0)+pow(y_vel_s,2.0)+pow(z_vel_s,2.0));
    	
    	/*find phase properties*/
    	den_l = C_R(cell, thread_l);/*the density of the liquid phase*/
    	vis_l = C_MU_EFF(cell, thread_l);/*the viscosity of the liquid phase*/
    	vol_s = C_VOF(cell, thread_s);/*solid volume fraction*/	
    	
    	/*particle diameter*/
    	diam_s = C_PHASE_DIAMETER(cell,thread_s);
    	
    	/*compute slip**/
    	abs_v = fabs(vel_l-vel_s);
    	
    	/*compute Reynolds number*/
    	rey = den_l*diam_s*abs_v/vis_l;
    	
    	/*compute drag coeff, dcef*/
    	if(rey <= 1000.0)
    	{
    		g0 = (1-vol_s)*rey;
    		g1 = 1.0+0.15*pow(g0,0.687);
    		dcef = 24.0*g1/g0;
    	}
    	else {dcef = 0.44;}
    	
    	/*compute fluid-solid exchange coefficient, k_l_s*/
    	g2 = pow((1.0-vol_s), -2.65);
    	k_wenyu = 0.75*dcef*(1-vol_s)*vol_s*den_l*abs_v*g2/diam_s;
    	g3 = (1-vol_s)*pow(diam_s, 2.0);
    	k_ergun = 150.0*vol_s*vol_s*vis_l/g3+1.75*vol_s*den_l*abs_v/diam_s;
    	
    	g4 = 262.5*(0.2-vol_s)/PAI;
    	stf = 0.5+atan(g4);/*stitching function*/
    	k_l_s = stf*k_ergun+(1.0-stf)*k_wenyu;
    	
    	return k_l_s;
    }
    

    不加载UDF,模型运算无误,加载UDF运算的时候会出现如下错误

    Updating solution at time level N... done.
      iter  continuity     u-water   u-hydrate     v-water   v-hydrate     w-water   w-hydrate     k-water   eps-water     time/iter
    # Divergence detected in AMG solver: mp-x-momentum -> Decreasing coarsening group size!
    # Divergence detected in AMG solver: mp-x-momentum -> Increasing relaxation sweeps!
    # Divergence detected in AMG solver: pressure correction -> Turning off correction scaling!
    # Divergence detected in AMG solver: pressure correction -> Increasing relaxation sweeps!
    
    Error: Divergence detected in AMG solver: pressure correction
    
    Error: Divergence detected in AMG solver: pressure correction
    Error Object: #f
    

    请教各位老师这是哪里出了问题。


  • Linux讲师 OpenFOAM讲师

    这就是发散了,这个 UDF return 的是k_l_s,这个值由几部分构成。
    这样,先 return 0 ,再慢慢往return里加项,看加到哪一项崩的。把那些项拆的越细越好,也可以在 return 之前加 printf 把各个项输出来,看看哪项的计算导致的发散,再去扣那项计算的具体过程。



  • 发散了。
    一般来说,都是你 return 的那一项太大了。可以先估算 return 的那一项的大小,然后按数量级缩小,直至不发散。



  • @bestucan 谢谢,按照您提供的思路找到了出错点



  • @zousiyu 找到的错误确实是return的那一项太大,非常感谢:146:


Log in to reply
 

CFD中文网 2016 - 2020 | 京ICP备15017992号-2