求助：一个两相流中曳力系数的UDF，可以interpret但计算出错

zoesjtuer

rt，下面是我在用户手册的例子的基础上改写的一个用于气液两相流的曳力系数udf，编译没有问题，但是开始计算50步左右就出现错误：Error: Floating point error: invalid number，请各位大侠看看udf哪里有问题呢？先谢谢大家了！
/*a programm for drag law*/
/*Cd=22.73*pow(Re,-0.849)*pow(Mo,0.020),0.5<=Re<=5*/
/*Cd=20.08*pow(Re,-0.636)*pow(Mo,0.046),5<=Re<=50*/
/*Mo=g*/
#include "udf.h"
#define Pi 1.*atan(1)
#define diam2 3.e-3
DEFINE_EXCHANGE_PROPERTY(drag_coef,cell,mix_thread,s_col,f_col)
{
Thread *thread_l, *thread_g;
real x_vel_g, x_vel_l, y_vel_g, y_vel_l, abs_v, slip_x, slip_y, density_g, density_l, viscosity_l, re, mo, sf_ts, drag_coef;
/*find the thread for the liquid(primary)and gas(secondary phase)*/
thread_l = THREAD_SUB_THREAD(mix_thread,s_col);/*liquid phase*/
thread_g = THREAD_SUB_THREAD(mix_thread,f_col);/*gas phase*/
/*find phase velocities and properties*/
x_vel_g = C_U(cell, thread_g);
y_vel_g = C_V(cell, thread_g);

x_vel_l = C_U(cell, thread_l);
y_vel_l = C_V(cell, thread_l);

slip_x = x_vel_l - x_vel_g;
slip_y = y_vel_l - y_vel_g;

density_l = C_R(cell,thread_l);
density_g = C_R(cell,thread_g);

viscosity_l = C_MU_L(cell, thread_l);

sf_ts = 0.00430;

/*compute slip*/

abs_v = sqrt(slip_x*slip_x + slip_y*slip_y);

/*compute re*/

re=density_l*abs_v*diam2/viscosity_l;

/*compute mo*/
mo = 9.81*pow(viscosity_l, 4)*(density_l - density_g)/(pow(sf_ts, 3)*pow(density_l, 2));

/*compute drag coefficient and return Cd*/
if(0.5<=re<=5)
   drag_coef = 22.73*pow(re, (-0.849))*pow(mo, 0.020);
else if(5<=re<=50)
   drag_coef = 20.08*pow(re, (-0.636))*pow(mo, 0.046);
return drag_coef;
}