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

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

[bestucan]

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

[zousiyu]

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

[推纱望月]

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

[推纱望月]

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

[彼岸]

可以把正确的UDF发来参考一下吗？我也做类似的复现研究，谢谢！