CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   Fluent UDF and Scheme Programming (http://www.cfd-online.com/Forums/fluent-udf/)
-   -   Help me to check my UDF (http://www.cfd-online.com/Forums/fluent-udf/117229-help-me-check-my-udf.html)

Liufeng_ustb May 5, 2013 05:48

Help me to check my UDF
 
hi, all.
i am eager to get your help.
it's my fist time to write my udf and i have many problems. followings are the main ones.
my simulation deals with heterogeneous reaction and i want to define the rate of heterogeneous reaction. u know, this can only be realized by udf for the fluent included in ansys 12.0.
for some reasons, i can use compiling udf, on the other hand, my udf s simple, so i chose inteprete it. however, when i finish all settings and start calculating, error comes out and shows "interpreted functions not supported". i try other fluent tutorials examples, they successfully execute. so i think my udf is wrong. here it is:

/*This UDF aims to define the reaction rate requied to simulate the reduction of pellet.
Here, gas is a primary phase mixture of three spieces: h2o and h2 and n2.
solid phase is also a mixture which consists of fe2o3 and its reduction product fe. */
/*本UDF是用来定义非均相反应速率的,用于模拟球团的还原过程;
这里的气相是一个由h2o,h2和n2组成的;
固相是由fe2o3和它的还原产物fe组成的*/
#include "udf.h"
#include "sg_mphase.h"
/*user inputs and some constants*/
/*用户输入的条件,其实就是在UDF关于Material and Phase中设置的情况,主要是顺序问题*/
#define MAX_SPE_EQNS_PRIM 3
/*total number of species in primary phase*/
/*primary phase中所含的所有物质的种类数目*/
#define index_reduction_primary 0
/*reduction species index in primary phase*/
/*为了在后面,满足UDF和求解器FLUENT交换数据的需要,需要对还原气体h2进行index的标记, 这个对应于在FLUENT设置mixture时的顺序问题*/
#define prim_index 0
/*index of primary phase*/
/*primary phase也需要用index来标记才能进行数据的交换*/
#define P_OPER 101325
/*operating pressure equal to GUI value*/
/*系统的操作压力是大气压,也就是说Bulk phase or primary phase or gas phase的总压,h2(g)的压力是其分压*/

/*end of user inputs*/
/*用户输入条件完毕,其实更应该理解为FLUENT中实现设定好的条件,尤其是对Mixture和Phas e的设置*/

/************************************************** ***********/
/* UDF for specifying an interfacial area density */
/*以下是本UDF真正的开始,本UDF对于are density
(应该是一个Cell中Secondary phase的particle的表面积)
的确定很关键,这是因为要用到表面积来描述反应的速率 */
/************************************************** ***********/
double p_equilibrium_h2(double t_solid,double p_gas)
/* */
/* Computes equilibrium pressure of hydrogen */
/* as function of temperature */
/* Equation is taken from THERMODYNAMIC THEORY IN SI, */
/* by Gibbs Equilibrium */
/* Returns pressure in PASCALS, given temperature in KELVIN */
/* */
/*氢气的平衡压里的计算要根据其与反应温度的关系来确定 */
/*公式是基于热力学的国际单位制 */
/*本质上是Gibbs用于化学平衡的等温方程 */
/*返回的压强单位是Pa,其中给出的温度使用的是热力学温标 */
/*根据反应动力学原理作出如下假设:
在fe2o3整个的逐级还原过程中,feo还原为金属fe这一步最困难,
在计算平衡常数和气相平衡浓度时,可以只考虑feo还原为fe的反应
feo+h2=fe+h2o */
{
double delt_g=10040-5.18*t_solid;
double interm=-delt_g/8.314/t_solid;
double p_equilibrium_h2=p_gas/(1+exp(interm));
return p_equilibrium_h2;
}


DEFINE_HET_RXN_RATE(diff_rate,c,t,hr,mw,yi,rr,rr_t )
{
Thread **pt=THREAD_SUB_THREADS(t);
/*pt is a pointer that stores all cell threads and every thread is also a pointer*/
/*This is a kind of way to get the pointers of subthread by using the known pointer of superthread*/
/*这里的pt是一个指针变量,就是一个数组,没一个元素也同时是以个指针,这个指针指向子线,每个单相用不 同的thread,
为了区分不同来自不同单相的thread,就必须用一个index来区分,这个就是子线指针所指向的内容, 就是*pt;
下面再总结一下,pt和*pt的含义:pt[i](其中i=0,1,2)是子线的指针,可以表示每个单相所用的thread;
*pt[i](其中i=0,1,2)是以个thread变量,表示的是cell组成的thread。*/
Thread *tp=pt[0];
Thread *ts=pt[1];
/*These two sentences corresponds primary phase and secondary phase
which are defined in FLUENT to them defined in UDF*/
/*这两条语句将前面所取得的subthread的指针分别给了两相,这样tp和ts就分别有了从FLUEN T中定义的两相的含义,
也可以这么理解,tp就是primary phase的thread,用来存储cell的信息。*/
int i;
real concentration_reduction_primary,accum = 0.,mole_frac_reduction_prim,
concentration_reduction_particle_surface;
/*Due to the rate of reduction which is related to the difference of h2 pressure between
in gas phase and at the surface of particle, we need to find the value of h2 pressure or concentration in
these two different places*/
/*由于还原反应的速率与h2在气相中的压力和在颗粒表面处的压力的差值有关,因此我们需要分别求出这两个地 方的压力,
压力也可以用h2的浓度来表示,在物理意义上是一致的*/
real T_prim=C_T(c,tp); /*primary phase (gas) temperature*/
real T_sec=C_T(c,ts); /*secondary phase (particle) temperature*/
/*UDF中关于Additional Macros for Writing UDFs in section 3有详细的介绍,对于Data Acess Macros这种类型的宏,
又有很多对我们有利的数据类型,例如Cell类型,C_T就可以告诉我们在某一个cell中的某一相的温度 */
real diam=C_PHASE_DIAMETER(c,ts);
/*secondary phase diameter*/
/*C_PHASE_DIAMETER和C_T一样,也是一个宏,可以得到在FLUENT中定义的某一相颗粒 的直径*/
real D_reduction_prim=3.96e-9;
/*primary phase effective species diffusivity*/
/*C_DIFF_EFF同样也是一个宏,传递了effective species diffusivity, D_reduction_prim 的单位是m2/s*/
real Re, Sc, Nu, urel, urelx,urely,urelz=0., mass_coeff, area_density,
flux_reduction ;
if(Data_Valid_P())
{
urelx = C_U(c,tp) - C_U(c,ts);
urely = C_V(c,tp) - C_V(c,ts);
#if RP_3D
urelz = C_W(c,tp) - C_W(c,ts);
#endif

urel = sqrt(urelx*urelx + urely*urely + urelz*urelz);
/*relative velocity*/
/*求相间的相对速度,是为了后面求出相关准数*/

Re = urel * diam * C_R(c,tp) / C_MU_L(c,tp);
/*雷诺数,就是用的我们所学的公式进行计算的,注意是primary ohase的密度和粘度,
需要注意的是,这里粘度用的是laminar viscosity并不是在求species of primary phase turbulent diffusivity,
我想可能是跟雷诺数的定义有关系吧*/

Sc = C_MU_L(c,tp) / C_R(c,tp) / D_reduction_prim ;
/*施密特数*/
Nu = 2. + 0.6 * pow(Re, 0.5)* pow(Sc, 0.333);
/*努赛尔数*/
mass_coeff = Nu * D_reduction_prim / diam ;
/*传质系数可以用努赛尔数得到,单位是m/s*/
for (i=0; i < MAX_SPE_EQNS_PRIM ; i++)
{
accum = accum + C_YI(c,tp,i)/mw[i][prim_index];
/*The function of variable accum is to acquire the total moles
in primary phase if given a 1kg primary phase*/
/*变量accum的含义:若给定1kg的primary phase,那么对于primary phase一共
就有accum moles的物质,引出accum就是为了把质量分数yi来换算成摩尔分数*/
}
mole_frac_reduction_prim = C_YI(c,tp,index_reduction_primary)
/ mw[index_reduction_primary][prim_index] / accum;
/*According the meaning of accum, we can readily know this formula can get the
mole fraction of reduction specie in primary phase if given 1kg primary phase*/
/*既然accum是1kg primary phase的摩尔数,那么1kg primary phase中h2的质量就可以
知道,其摩尔数也就可以求得。因此水的摩尔分数自然也就很容易得到了*/
concentration_reduction_primary = mole_frac_reduction_prim * C_P(c,tp)
/ 8.314 / T_prim ;
/*This formula is an application of P=cRT, besides, Dalton's law is also used for
p(h2) in primary phase*/
/*该式应用了理想气体的状态方程求出了h2在气相中的浓度,单位是mol/m3,在求h2分压的
时候,还使用的道尔顿分压定律*/
concentration_reduction_particle_surface = p_equilibrium_h2(T_sec,C_P(c,tp))/8.314/T_sec ;
/*This formula is much easier because at the surface of particle, the h2 pressure is
exactly the equlibrium pressure, but we need to notice the temperature is the temperature of
secondary phase*/
/*这个式子还是理想气体状态方程,只不过这次的压力就是颗粒表面的平衡压力,此外,温度也使用的
是secondary phase的温度*/
area_density = 6. * C_VOF(c,ts) / diam ;
/*area_density is a variable representing the area of the particle of secondary phase in a cell
which also includes primary phase. It will be used to get rr*/
/*area_density得到了单位体积,也就是一个cell中secondary phase的表面积,其单位是/m*/
flux_reduction = mass_coeff *
(concentration_reduction_particle_surface - concentration_reduction_primary ) ;
/*flux_reduction is rate of trasporting mass, its theoretical base is obvious for us*/
/*flux_reduction就是单位面积的传质速率,对于这个还原过程来说,就是我们简化的反应的速率 ,但是其单位是mol/(m2*s)*/
*rr = area_density * flux_reduction ;
/* *rr is a pointer, its entry is exactly the reaction rate for hetergeneous reaction*/
/* *rr就是我们要得到的反应速率,单位是mol/(m3*s)*/
}
}

vasava May 7, 2013 07:32

I think you are not applying the reaction rate to the whole domain (I assume thats what you want to do). Instead you are using 'Thread' and limiting it to a face or a boundary. I suggest you consider using 'Domain' instead.

Please read the help section before you implement my suggestion.

Liufeng_ustb May 7, 2013 11:25

Thanks for your advice. I will have a try tomorrow using the computer in experiment room. If I have another trouble, I will turn you for help.
Thank you!


All times are GMT -4. The time now is 09:03.