本案例演示利用Fluent UDF自定义欧拉多相流模型中的颗粒阻力律以模拟流化床中的流体流动。

1 问题描述

在欧拉多相流模型中常使用Syamlal-O’Brien模型描述颗粒与连续相之间的相互作用力。

Fluent中默认的流体-固体颗粒阻力系数描述为：

其中为固相颗粒的体积分数；为固相颗粒的密度；为固相颗粒的弛豫时间，其定义为：

其中，为连续相的动力粘度；为颗粒相的粒径。

式(1)中的用于区分不同的模型，对于Syamlal_O’Brien模型，函数定义为：

其中为连续相的体积分数；阻力系数定义为：

式中为颗粒雷诺数，其定义为：

式中为固相的末端速度系数。

式中，；当，当。

在需如下图所示的流化床模拟中，利用默认Syamlal-O’Brien模型的默认常数0.8及2.65预测得到的最小流态化速度为21 cm/s，而实验观察得到的最小流态化速度为8 cm/s，因此需要对默认的Syamlal-O'Brien模型中的常数进行修改，经过调整后常数分别修改为0.281632及9.07696。

得到末端速度系数后，即可利用上面的公式回代得到系数。

计算模型如下图所示，计算区域尺寸为1 m ×0.15 m，空气以0.25 m/s速度进入到计算区域，顶部边界为出口。流化床初始条件下填充率为0.55。

2 准备UDF

根据前面的描述，编写UDF如下所示。

#include "udf.h"
#include "sg_mphase.h"
 
#define pi 4. * atan(1.)
#define diam2 3.e-4
 
DEFINE_EXCHANGE_PROPERTY(custom_drag_syam, cell, mix_thread, s_col, f_col)
{
     Thread *thread_g, *thread_s;
     real x_vel_g, x_vel_s, y_vel_g, y_vel_s, abs_v, slip_x, slip_y,
         rho_g, rho_s, mu_g, reyp, afac,
         bfac, void_g, vfac, fdrgs, taup, k_g_s;
 
     /* 根据主相和次相返回各自的Thread
        s_col返回的是气相，f_col返回颗粒相*/
     thread_g = THREAD_SUB_THREAD(mix_thread, s_col);
     thread_s = THREAD_SUB_THREAD(mix_thread, f_col);
 
     /*得到各相的速度与物性*/
     x_vel_g = C_U(cell, thread_g);
     y_vel_g = C_V(cell, thread_g);
     x_vel_s = C_U(cell, thread_s);
     y_vel_s = C_V(cell, thread_s);
     slip_x = x_vel_g - x_vel_s;
     slip_y = y_vel_g - y_vel_s;
     rho_g = C_R(cell, thread_g);
     rho_s = C_R(cell, thread_s);
     mu_g = C_MU_L(cell, thread_g);
 
     /*计算得到滑移速度的合速度*/
     abs_v = sqrt(slip_x * slip_x + slip_y * slip_y);
     /*得到颗粒雷诺数*/
     reyp = rho_g * abs_v * diam2 / mu_g;
     /*计算得到颗粒弛豫时间 */
     taup = rho_s * diam2 * diam2 / 18. / mu_g;
     /*得到颗粒体积分数*/
     void_g = C_VOF(cell, thread_g);
 
     /*得到系数bfac*/
     afac = pow(void_g, 4.14);
     if (void_g <= 0.85)
          bfac = 0.281632 * pow(void_g, 1.28);
     else
          bfac = pow(void_g, 9.076960);
    /*代公式计算*/
     vfac = 0.5 * (afac - 0.06 * reyp + sqrt(0.0036 * reyp * reyp + 0.12 * reyp * (2. * bfac - afac) + afac * afac));
     fdrgs = void_g * (pow((0.63 * sqrt(reyp) / vfac + 4.8 * sqrt(vfac) / vfac), 2)) / 24.0;
     k_g_s = (1. - void_g) * rho_s * fdrgs / taup;
 
     return k_g_s;
}