气液两相流源项UDF

iwangchunyu

在处理气液两相流源项的问题上，我采用了下面的方程(见本楼底附件)：

在编写程序的时候，有两个问题：
1. 体积分数梯度和温度梯度的乘积怎么表示？
2. ds[eqn]项如何编写？

目前我是这样做的：
1. 两个梯度的乘积，表示为向量的点积，即各自的分量乘积和；
2. 主相气相中的ds[eqn]为对方程中气相分数afa_v求导，气相分数的梯度看做常数（这里感觉有问题），温度梯度、导热系数都是常数；
3. 第二相液相中ds[eqn]为对方程中液相分数afa_l求导，同样，气相分数的梯度、温度梯度哦、导热系数也为常数。

按照这样的思路，编写了一下udf，读入成功，初始化或计算时报错（ACCESS_VIOLATION)。考虑问题的原因，可能是第一步计算前，梯度无法获得。在未导入udf的情况下，计算数步，再导入udf计算，依然报错。看来还是程序的问题。
希望论坛里的各位高手可以指点一二。

udf如下：
#include "udf.h"  

#include "sg_mphase.h"  
#define T_SAT 319   
#define LAT_HT 2.392e6  

DEFINE_SOURCE(vap_src, cell, pri_th, dS, eqn)
{
Thread * mix_th, *sec_th;
real m_dot_v,q;
mix_th = THREAD_SUPER_THREAD(pri_th);
sec_th = THREAD_SUB_THREAD(mix_th, 1);
if(C_T(cell, mix_th)<=T_SAT)
{
q=-((C_K_L(cell,sec_th)*C_VOF(cell,sec_th)*(C_T_G(cell,sec_th)[0]*C_VOF_G(cell,pri_th)[0]+C_T_G(cell,sec_th)[1]*C_VOF_G(cell,pri_th)[1]))-(C_K_L(cell,pri_th)*C_VOF(cell,pri_th)*(C_T_G(cell,pri_th)[0]*C_VOF_G(cell,pri_th)[0]+C_T_G(cell,pri_th)[1]*C_VOF_G(cell,pri_th)[1])));
m_dot_v=q/LAT_HT;  
dS[eqn] = -C_K_L(cell,pri_th)*(C_T_G(cell,pri_th)[0]*C_VOF_G(cell,pri_th)[0]+C_T_G(cell,pri_th)[1]*C_VOF_G(cell,pri_th)[1])/LAT_HT;
}
else
{
q=((C_K_L(cell,sec_th)*C_VOF(cell,sec_th)*(C_T_G(cell,sec_th)[0]*C_VOF_G(cell,pri_th)[0]+C_T_G(cell,sec_th)[1]*C_VOF_G(cell,pri_th)[1]))-(C_K_L(cell,pri_th)*C_VOF(cell,pri_th)*(C_T_G(cell,pri_th)[0]*C_VOF_G(cell,pri_th)[0]+C_T_G(cell,pri_th)[1]*C_VOF_G(cell,pri_th)[1])));
m_dot_v=q/LAT_HT;
dS[eqn] =dS[eqn] = C_K_L(cell,pri_th)*(C_T_G(cell,pri_th)[0]*C_VOF_G(cell,pri_th)[0]+C_T_G(cell,pri_th)[1]*C_VOF_G(cell,pri_th)[1])/LAT_HT;
}
return m_dot_v;
}

DEFINE_SOURCE(liq_src, cell, sec_th, dS, eqn)
{
Thread *mix_th, *pri_th;  
  real m_dot_l,  q;
  mix_th = THREAD_SUPER_THREAD(sec_th);  
  pri_th = THREAD_SUB_THREAD(mix_th,0);
  if(C_T(cell, mix_th)<=T_SAT)  
{
q=((C_K_L(cell,sec_th)*C_VOF(cell,sec_th)*(C_T_G(cell,sec_th)[0]*C_VOF_G(cell,pri_th)[0]+C_T_G(cell,sec_th)[1]*C_VOF_G(cell,pri_th)[1]))-(C_K_L(cell,pri_th)*C_VOF(cell,pri_th)*(C_T_G(cell,pri_th)[0]*C_VOF_G(cell,pri_th)[0]+C_T_G(cell,pri_th)[1]*C_VOF_G(cell,pri_th)[1])));
m_dot_l=q/LAT_HT;  
dS[eqn] =C_K_L(cell,sec_th)*(C_T_G(cell,sec_th)[0]*C_VOF_G(cell,pri_th)[0]+C_T_G(cell,sec_th)[1]*C_VOF_G(cell,pri_th)[1])/LAT_HT;

}
else
{
q=-((C_K_L(cell,sec_th)*C_VOF(cell,sec_th)*(C_T_G(cell,sec_th)[0]*C_VOF_G(cell,pri_th)[0]+C_T_G(cell,sec_th)[1]*C_VOF_G(cell,pri_th)[1]))-(C_K_L(cell,pri_th)*C_VOF(cell,pri_th)*(C_T_G(cell,pri_th)[0]*C_VOF_G(cell,pri_th)[0]+C_T_G(cell,pri_th)[1]*C_VOF_G(cell,pri_th)[1])));
m_dot_l=q/LAT_HT;  
dS[eqn] =-C_K_L(cell,sec_th)*(C_T_G(cell,sec_th)[0]*C_VOF_G(cell,pri_th)[0]+C_T_G(cell,sec_th)[1]*C_VOF_G(cell,pri_th)[1])/LAT_HT;
}
return m_dot_l;
}

DEFINE_SOURCE(enrg_src, cell, mix_th, dS, eqn)
{
  Thread *pri_th, *sec_th;
  real m_dot, q;
  pri_th = THREAD_SUB_THREAD(mix_th, 0);
  sec_th = THREAD_SUB_THREAD(mix_th, 1);
  if(C_T(cell, mix_th)<=T_SAT)
{
q=((C_K_L(cell,sec_th)*C_VOF(cell,sec_th)*(C_T_G(cell,sec_th)[0]*C_VOF_G(cell,pri_th)[0]+C_T_G(cell,sec_th)[1]*C_VOF_G(cell,pri_th)[1]))-(C_K_L(cell,pri_th)*C_VOF(cell,pri_th)*(C_T_G(cell,pri_th)[0]*C_VOF_G(cell,pri_th)[0]+C_T_G(cell,pri_th)[1]*C_VOF_G(cell,pri_th)[1])));
dS[eqn] =0;
}
else
{
q=-((C_K_L(cell,sec_th)*C_VOF(cell,sec_th)*(C_T_G(cell,sec_th)[0]*C_VOF_G(cell,pri_th)[0]+C_T_G(cell,sec_th)[1]*C_VOF_G(cell,pri_th)[1]))-(C_K_L(cell,pri_th)*C_VOF(cell,pri_th)*(C_T_G(cell,pri_th)[0]*C_VOF_G(cell,pri_th)[0]+C_T_G(cell,pri_th)[1]*C_VOF_G(cell,pri_th)[1])));
dS[eqn] =0;
}
return q;  
}

[ 本帖最后由 iwangchunyu 于 2013-12-24 15:23 编辑 ]

12535111867

帮楼主顶下，我是个菜虫，希望楼主做好，分享一下经验