Euler双流体模型修改曳力模型后发散

[zousiyu]

fluent默认模型收敛，修改后直接发散。提示：

Updating solution at time level N... done.
  iter  continuity     u-water       u-air     v-water       v-air      vf-air     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!
# You may try the enhanced divergence recovery with (rpsetvar 'amg/protective-enhanced? #t)
# Divergence detected in AMG solver: pressure correction -> Turning off correction scaling!
# Divergence detected in AMG solver: pressure correction -> Increasing relaxation sweeps!
# You may try the enhanced divergence recovery with (rpsetvar 'amg/protective-enhanced? #t)

Error: Divergence detected in AMG solver: pressure correction

Error: Divergence detected in AMG solver: pressure correction
Error Object: #f

这种情况是因为曳力突然为0了，或者无穷大？？？？？那位老司机教教我

[李东岳]

怎么修改的模型？

[zousiyu]

@东岳 计算了新的曳力系数

#include "udf.h"
 
DEFINE_EXCHANGE_PROPERTY(F_gl, c, mixture_thread, second_col, first_col)
{
    /* coefficient and parameters */
    real E_1 = 180;
    real E_2 = 1.8;
    real r_p = 0.0015;
    real d_p = 2 * r_p;
 
    /* declaration */
    Thread *thread_g, *thread_l;
    real x_vel_g, y_vel_g;
    real x_vel_l, y_vel_l;
    real slip_x, slip_y, abs_slip_vel;
    real rho_g, rho_l, mu_g, mu_l;
    real void_g, void_l, void_s;
    real K_gl, tmp_1, tmp_2;
 
    /* find the threads for the gas and liquid */
    thread_g = THREAD_SUB_THREAD(mixture_thread, second_col); /* gas */
    thread_l = THREAD_SUB_THREAD(mixture_thread, first_col);  /* liquid */
 
    /* get phase velocity and fraction */
    x_vel_g = C_U(c, thread_g);
    y_vel_g = C_V(c, thread_g);
    x_vel_l = C_U(c, thread_l);
    y_vel_l = C_V(c, thread_l);
    slip_x = x_vel_g - x_vel_l;
    slip_y = y_vel_g - y_vel_l;
    mu_g = C_MU_L(c, thread_g);
    mu_l = C_MU_L(c, thread_l);
    rho_g = C_R(c, thread_g);
    rho_l = C_R(c, thread_l);
    void_g = C_VOF(c, thread_g);
    void_l = C_VOF(c, thread_l);
    void_s = 0.365;
 
    /* slip velocity */
    abs_slip_vel = sqrt(slip_x * slip_x + slip_y * slip_y);
 
    /* compute drag and return drag coeff */
    K_gl = (void_g / (void_g + void_l)) * (tmp_1 + tmp_2);
    tmp_1 = (E_1 * mu_g * pow(1 - void_g, 2) / pow(void_g, 2) / pow(d_p, 2)) *
            pow((void_s / (1 - void_g)), 3 / 2);
    tmp_2 = (E_2 * rho_g * abs_slip_vel * (1 - void_g) / void_g / d_p) *
            pow((void_s / (1 - void_g)), 1 / 3);
    return K_gl;
}

[zousiyu]

@东岳 老师有什么意见吗？

[李东岳]

不懂Fluent，不过43行代码之前

K_gl = (void_g / (void_g + void_l)) * (tmp_1 + tmp_2);

是不是放后面比较好？然后void_g + void_l不应该等于1么？

[zousiyu]

@东岳 谢谢，这个错误我已经发现了。问题已经解决了，因为void_g会为1，那么1-void_g就是0，然后会产生NaN错误。
我直接把void_g为0的情况强制修改为0.5了

[李东岳]

我直接把void_g为0的情况强制修改为0.5了

我更喜欢变成 void_g = void_g + 1e-6

[walten]

@东岳 同意，我一般都是改成 + 1e-08