找回密码
 注册
查看: 3299|回复: 3

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

[复制链接]
发表于 2012-9-6 14:08:14 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?注册

x
各位好!
    最近在编写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 编辑 ]
帮助文件面宏说明.jpg
表3.2.20.jpg
表3.2.21.jpg
表3.2.22.jpg
发表于 2012-9-6 16:56: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_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的平均值来求。
 楼主| 发表于 2012-9-6 20:20:57 | 显示全部楼层
原帖由 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 编辑 ]
发表于 2014-10-6 16:59:08 | 显示全部楼层

回复 1# JohnQiang 的帖子

我想你是不是可以用一下C_Perfix呢?
您需要登录后才可以回帖 登录 | 注册

本版积分规则

快速回复 返回顶部 返回列表