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



  • 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了,或者无穷大?????那位老司机教教我

    0_1541423078443_1d6ba055-4e99-4c56-8500-269256ebe3c7-image.png


  • 网格教授 OpenFOAM教授 管理员

    怎么修改的模型?



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

    #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;
    }
    


  • @东岳 老师有什么意见吗?


  • 网格教授 OpenFOAM教授 管理员

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

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

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



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


  • 网格教授 OpenFOAM教授 管理员

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

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