求助，udf编写yplus

kxflove1991

我想通过UDF自定义一个湍流普朗特特数，具体的函数关系如图所示：

但是无法复现文献中的结果，现在怀疑yplus的编写有问题，采用UDF计算时非常容易发散，请各位专家给点建议。代码如下所示：
计算区域为二位轴对称的圆管，半径为4mm。
/* Specifying a Variable Temperature Prandtl Number */
#include "udf.h"
#include "math.h"
DEFINE_PRANDTL_T(Prt,c,t)
{
   real pr_t;                                                           /* turbulent Prandtl */
    real mu = C_MU_L(c,t);                                   /* molecular dynamic viscosity */
    real mu_t = C_MU_T(c,t);                                 /* turbulent dynamic viscosity */
    real thermalCon=C_K_L(c,t);                              /* molecular thermal conductivity */
    real cp=C_CP(c,t);                                               /* specific heat */
    real pr_mol=mu*cp/thermalCon;                    /* molecular pr */
    real density=C_R(c,t);                                               /* molecular density */
    real mu_ratio=mu_t/mu;                                   /* mu_ratio */
    real velocity_dudy=C_DUDY(c,t);                          /* velocity derivative */
    real Tao_tur=mu*velocity_dudy;                           /* turbulence shear stress */
    real velocity_F=sqrt(fabs(Tao_tur/density));  /* Friction Velocity */
    real y;
    real y_plus;
    real x[2];                                      
    C_CENTROID(x,c,t);
    y=x[1];               
   y_plus=density*velocity_F*(0.004-y)/mu;
               
                if  (y_plus<=5)
                        pr_t=1.0;
                else if (y_plus<100 && y_plus>5)
                        pr_t=0.5+0.1*pr_mol/mu_ratio;
                else
                        pr_t=0.85;       
               
        return pr_t;
}