请教：F_U(f,t)等面宏不能用在密度基求解器上

JohnQiang

各位好！
    最近在编写UDF程序，编译连接都没问题，可是计算中遇到了问题(这个问题耗了一个多月，还没解决)：计算到UDS方程的对流项时出错。检查分析发现是宏F_U(f,t)、F_V(f,t)和F_W（f,t）的问题，查找原因发现宏F_U(f,t)、F_U(f,t)和F_W（f,t）只能用在压力基求解器（Pressure Based Solver）情况下（附件为UDF帮助文件的说明）。
    关键是我的计算求解器必须用密度基求解器（Density Based Solver），因为用到了Real-Gas-Model。
    这种情况该如何处理呢？
    还望各位不吝赐教，谢谢！

下面是我的程序段：（程序不是重点，只是想说明是用Message（...）检查出问题的）
DEFINE_UDS_FLUX(x_motion_flux,f,t,i)
{
cell_t c0;
cell_t c1=-1;
Thread *t0;
Thread *t1=NULL;
real NV_VEC(psi_vec);
real NV_VEC(A);
real flux=0.0;
c0=F_C0(f,t);
t0=F_C0_THREAD(f,t);
F_AREA(A,f,t);
if(BOUNDARY_FACE_THREAD_P(t)) /* Most face values will be available*/
{
        real dens;
Message("\n x-momentum-flux     001  \n");
        /* Depending on its BC, density may not be set on face thread. */

        if(NNULLP(THREAD_STORAGE(t,SV_DENSITY)))
                dens=F_R(f,t)*C_UDSI(c0,t0,b_continuty); /* Set dens to face value if available, */
        else
                dens=C_R(c0,t0)*C_UDSI(c0,t0,b_continuty); /* else set dens to cell value */
Message("\nx-momentum-flux  F-density=%e \n",dens);
Message("\n22222222222222222\n");
Message("\n f-u=%e f-v=%e f-w=%e \n",F_U(f,t),F_V(f,t),F_W(f,t));
Message("\n333333333333333333333333\n");       
NV_DS(psi_vec,=,F_U(f,t),F_V(f,t),F_W(f,t),*,dens);
        flux=NV_DOT(psi_vec,A); /* Flux through Face. */
Message("\n\n\n\n FLux=%e\n",flux);
}
else
{
        c1=F_C1(f,t); /* Get cell on other side of face. */
        t1=F_C1_THREAD(f,t);
Message("\n x-momentum-flux     002  \n");        NV_DS(psi_vec,=,C_U(c0,t0),C_V(c0,t0),C_W(c0,t0),*,C_R(c0,t0)*C_UDSI(c0,t0,b_continuty));
        NV_DS(psi_vec,+=,C_U(c1,t1),C_V(c1,t1),C_W(c1,t1),*,C_R(c1,t1)*C_UDSI(c1,t1,b_continuty));
        flux=NV_DOT(psi_vec,A)/2.0; /* Average flux through face. */
}
Message("\n x_motion_flux=%e \n",flux);
    return flux;
}

[ 本帖最后由 JohnQiang 于 2012-9-6 14:32 编辑 ]

gearboy78

举个例子，如果算x方向速度。最简单的就是直接用c0和c1值的平均值
face_vel=0.5*(C_U(c0,t0)+C_U(c1,t1));
如果要精度高一些，就用一阶泰勒展开的平均值
face_vel0=C_U(c0,t0)+NV_DOT(C_U_G(c0,t0),dr0);
face_vel1=C_U(c1,t1)+NV_DOT(C_U_G(c1,t1),dr1);
face_vel=0.5*(face_vel0+face_vel1);
其中dr0是连接c0到face的矢量，dr1是连接c1到face的矢量.可用INTERIOR_FACE_GEOMETRY获取。

如果是边界面，由于不存在c1，就直接用c0的值或者c0的一阶泰勒展开。

此外也可用面的node的平均值来求。

JohnQiang

原帖由 gearboy78 于 2012-9-6 16:56 发表

                               
登录/注册后可看大图

举个例子，如果算x方向速度。最简单的就是直接用c0和c1值的平均值
face_vel=0.5*(C_U(c0,t0)+C_U(c1,t1));
如果要精度高一些，就用一阶泰勒展开的平均值
face_vel0=C_U(c0,t0)+NV_DOT(C_U_G(c0,t0),dr0);
face_v ...

谢谢gearboy78的指点，有点豁然开朗的感觉了。
感觉是这样的，还是自己数学或者计算流体力学没学好，用理论解决问题的能力差，没往深处想，我还是应该好好琢磨一下它程序的根本的...学以致用，我必须得反思...
总之，谢谢gearboy78

[ 本帖最后由 JohnQiang 于 2012-9-6 20:26 编辑 ]

sam_paris

我想你是不是可以用一下C_Perfix呢？