小弟最近在做电弧双温模型，两个能量方程都是通过uds写的，程序写好后导入fluent计算几千步就发散了，速度入口附近的速度和压力变得非常大，怀疑是对流项的问题，但是不太会修改，请各位大大帮忙看一下，万分感激！小木虫里不会发图片，发不了公式，这里先附上我的代码：
DEFINE_UDS_FLUX(uds3_flux, f, t, i)
{
        cell_t c0, c1 = -1;
        Thread *t0, *t1 = NULL;
        real NV_VEC(psi_vec);
        real NV_VEC(V);
        real NV_VEC(A), flux = 0.0;
        real ke = 1.38e-23;
        c0 = F_C0(f, t);
        t0 = F_C0_THREAD(f, t);        //内部面，气体流入的为C0，流出为C1。边界面上只有C0没有C1
        F_AREA(A, f, t);
        /* If face lies at domain boundary, use face values; */
        /* If face lies IN the domain, use average of adjacent cells. */
        if (BOUNDARY_FACE_THREAD_P(t)) /*Most face values will be available*/
        {
                real dens;

                /* Depending on its BC, density may not be set on face thread*/
                if (NNULLP(THREAD_STORAGE(t, SV_DENSITY)))
                {
                        dens = edensity(F_UDSI(c0, t0, 3), F_UDSI(c0, t0, 4));  /* Set dens to face value if available */
                        C_UDMI(c0, t0, 16) = edensity(F_UDSI(c0, t0, 3), F_UDSI(c0, t0, 4));
                }
                else
                {
                        dens = edensity(C_UDSI(c0, t0, 3), C_UDSI(c0, t0, 4)); /* else, set dens to cell value */
                        C_UDMI(c0, t0, 17) = edensity(C_UDSI(c0, t0, 3), C_UDSI(c0, t0, 4));
                }
                psi_vec[0] = F_U(f, t)*2.5*dens*ke;
                psi_vec[1] = F_V(f, t)*2.5*dens*ke;
                /*NV_DS(psi_vec, =, F_U(f,t), F_V(f,t), F_W(f,t), *, dens*ke);*/
                flux = NV_DOT(psi_vec, A); /* flux through Face */
               
        }
        else
        {

                c1 = F_C1(f, t);   /* Get cell on other side of face */
                t1 = F_C1_THREAD(f, t);
                psi_vec[0] = C_U(c0, t0)*2.5*ke*edensity(C_UDSI(c0, t0, 3), C_UDSI(c0, t0, 4));
               
                C_UDMI(c0, t0, 18) = C_U(c0, t0)*2.5*ke*edensity(C_UDSI(c0, t0, 3), C_UDSI(c0, t0, 4));
                C_UDMI(c1, t1, 19) = C_V(c1, t1)*2.5*ke*edensity(C_UDSI(c1, t1, 3), C_UDSI(c1, t1, 4));
                psi_vec[1] = C_V(c1, t1)*2.5*ke*edensity(C_UDSI(c1, t1, 3), C_UDSI(c1, t1, 4));
                flux = NV_DOT(psi_vec, A)/2.0; /* Average flux through face */
                C_UDMI(c0, t0, 23) = flux;
               
        }
        /* ANSYS Fluent will multiply the returned value by phi_f (the scalar’s
        value at the face) to get the ‘‘complete’’ advective term. */
        C_UDMI(c0, t0, 24) = flux;
        return flux;
}@wuming524 返回小木虫查看更多