求助两相流udf问题

tianhonghh

请各位大神帮我看看这个udf，计算不了几步就发散
#include "udf.h"
#include "sg_mphase.h"
DEFINE_SOURCE(liq_src,c,pri_th,dS,eqn)
{
Thread *mix_th,*sec_th;
real m_dot_l;
real stc;
real fv;
real pres_vap;
real pres_sat;
real pres_turb;
mix_th=THREAD_SUPER_THREAD(pri_th);
sec_th=THREAD_SUB_THREAD(mix_th,1);
fv=C_VOF(c,sec_th)*C_R(c,sec_th)/(C_VOF(c,sec_th)*C_R(c,sec_th)+C_VOF(c,pri_th)*C_R(c,pri_th));
stc=5.98633e-7*C_T(c,mix_th)*C_T(c,mix_th)-3.21398e-4*C_T(c,mix_th)+0.03015;
pres_sat=0.04782*C_T(c,mix_th)*C_T(c,mix_th)*C_T(c,mix_th)*C_T(c,mix_th)-4.33987*C_T(c,mix_th)*C_T(c,mix_th)*C_T(c,mix_th)-172.30441*C_T(c,mix_th)*C_T(c,mix_th)+27971.53714*C_T(c,mix_th)-734733.4355;
pres_turb=0.39*C_R(c,mix_th)*C_K(c,mix_th);
pres_vap=pres_sat+0.5*pres_turb;
if(C_P(c,mix_th)<=pres_vap)
{
m_dot_l=-0.02*sqrt(C_K(c,mix_th))/stc*C_R(c,pri_th)*C_R(c,sec_th)*sqrt((2*(pres_vap-C_P(c,mix_th)))/(3*C_R(c,pri_th)))*(1-fv);
dS[eqn]=-0.02*sqrt(C_K(c,mix_th))/stc*C_R(c,pri_th)*C_R(c,sec_th)*sqrt((2*(pres_vap-C_P(c,mix_th)))/(3*C_R(c,pri_th)));
}
else
{
m_dot_l=0.01*sqrt(C_K(c,mix_th))/stc*C_R(c,pri_th)*C_R(c,pri_th)*sqrt((2*(C_P(c,mix_th)-pres_vap))/(3*C_R(c,pri_th)))*fv;
dS[eqn]=0;
}
return m_dot_l;
}

DEFINE_SOURCE(vap_src,c,sec_th,dS,eqn)
{
Thread *mix_th,*pri_th;
real m_dot_v;
real stc;
real fv;
real pres_vap;
real pres_sat;
real pres_turb;
mix_th=THREAD_SUPER_THREAD(sec_th);
pri_th=THREAD_SUB_THREAD(mix_th,0);
fv=C_VOF(c,sec_th)*C_R(c,sec_th)/(C_VOF(c,sec_th)*C_R(c,sec_th)+C_VOF(c,pri_th)*C_R(c,pri_th));
stc=5.98633e-7*C_T(c,mix_th)*C_T(c,mix_th)-3.21398e-4*C_T(c,mix_th)+0.03015;
pres_sat=0.04782*C_T(c,mix_th)*C_T(c,mix_th)*C_T(c,mix_th)*C_T(c,mix_th)-4.33987*C_T(c,mix_th)*C_T(c,mix_th)*C_T(c,mix_th)-172.30441*C_T(c,mix_th)*C_T(c,mix_th)+27971.53714*C_T(c,mix_th)-734733.4355;
pres_turb=0.39*C_R(c,mix_th)*C_K(c,mix_th);
pres_vap=pres_sat+0.5*pres_turb;
if(C_P(c,mix_th)<=pres_vap)
{
m_dot_v=0.02*sqrt(C_K(c,mix_th))/stc*C_R(c,pri_th)*C_R(c,sec_th)*sqrt((2*(pres_vap-C_P(c,mix_th)))/(3*C_R(c,pri_th)))*(1-fv);
dS[eqn]=0;
}
else
{
m_dot_v=-0.01*sqrt(C_K(c,mix_th))/stc*C_R(c,pri_th)*C_R(c,pri_th)*sqrt((2*(C_P(c,mix_th)-pres_vap))/(3*C_R(c,pri_th)))*fv;
dS[eqn]=-0.01*sqrt(C_K(c,mix_th))/stc*C_R(c,pri_th)*C_R(c,pri_th)*sqrt((2*(C_P(c,mix_th)-pres_vap))/(3*C_R(c,pri_th)));
}
return m_dot_v;
}

DEFINE_SOURCE(enrg_src,c,mix_th,dS,eqn)
{
Thread *pri_th,*sec_th;
real m_dot;
real hc;
real source;
real stc;
real fv;
real pres_vap;
real pres_sat;
real pres_turb;
pri_th=THREAD_SUB_THREAD(mix_th,0);
sec_th=THREAD_SUB_THREAD(mix_th,1);
fv=C_VOF(c,sec_th)*C_R(c,sec_th)/(C_VOF(c,sec_th)*C_R(c,sec_th)+C_VOF(c,pri_th)*C_R(c,pri_th));
stc=5.98633e-7*C_T(c,mix_th)*C_T(c,mix_th)-3.21398e-4*C_T(c,mix_th)+0.03015;
pres_sat=0.04782*C_T(c,mix_th)*C_T(c,mix_th)*C_T(c,mix_th)*C_T(c,mix_th)-4.33987*C_T(c,mix_th)*C_T(c,mix_th)*C_T(c,mix_th)-172.30441*C_T(c,mix_th)*C_T(c,mix_th)+27971.53714*C_T(c,mix_th)-734733.4355;
pres_turb=0.39*C_R(c,mix_th)*C_K(c,mix_th);
pres_vap=pres_sat+0.5*pres_turb;
hc=-3.01379e-4*C_T(c,mix_th)*C_T(c,mix_th)*C_T(c,mix_th)+0.05766*C_T(c,mix_th)*C_T(c,mix_th)-4.80036*C_T(c,mix_th)+364.96837;
if(C_P(c,mix_th)<=pres_vap)  
{
m_dot=-0.02*sqrt(C_K(c,mix_th))/stc*C_R(c,pri_th)*C_R(c,sec_th)*sqrt((2*(pres_vap-C_P(c,mix_th)))/(3*C_R(c,pri_th)))*(1-fv);
source=m_dot*hc;
dS[eqn]=0;
}
else
{
m_dot=0.01*sqrt(C_K(c,mix_th))/stc*C_R(c,pri_th)*C_R(c,pri_th)*sqrt((2*(C_P(c,mix_th)-pres_vap))/(3*C_R(c,pri_th)))*fv;
source=m_dot*hc;
dS[eqn]=0;
}
return source;
}

langlovezhuzhu

请问你的问题解决了吗？我的也有类似问题、、、

tianhonghh

还没有

gearboy

为啥满世界都是这个烂模型，试想如下程序
if(a<=0) a=1
if(a>=1) a=0;
那a的值就会在0和1之间来回蹦，刚要到0又被你改成1，刚到1又被你改成0，残差还能稳定下降吗？

tianhonghh

您好。请问该怎么改呢？我的判定条件是压力，根据压力的关系确定质量表达式，即
if(p<pv) m=....;if(p>pv) m=....跟楼主说的if(a<=0) a=1 if(a>=1) a=0一样吗？

gearboy

差不多类似吧，试想p>pv时，假设m=1，p<pv时m=-1，那么当p>pv时你的源项就是m=1，必然引起系统变化，残差上升，然后一步一步迭代，好不容易下降稳定下来，这时候又变成p<pv了，系统必然又不稳定，残差再次上升。然后又开始迭代逐渐稳定，刚稳定又出现p>pv。所以整个过程残差会上下跳。除非在p=pv的地方，你的两个m的计算结果是相同的，也就是满足连续性。可惜，我看你的两个m的表达式计算结果是不同的。

[ 本帖最后由 gearboy 于 2014-4-3 06:42 编辑 ]

sanhuzb

楼主，你的问题解决了吗？我现在也碰到这个问题，求解答，谢谢

tianhonghh

没有解决呢