[m.r.soufivand]

As in the article " https://doi.org/10.1016/j.cep.2019.04.015 ", I want to simulate direct condensation of steam to subcooled water spray.

The results are good when I simulate only hydrodynamics, but when I compile the UDF to mass transfer and add source terms, In the first few iterations, it diverges and gives the following error:
"WARNING: Invalid cp (0.000000e+00 J/kgK) for water-vapor at temperature nan K
Error at host: floating point exception
Error: floating point exception
Error Object: #f"



I have attached the UDF to this post. I would be grateful if you could help me.

[AlexanderZ]

your code

#include "udf.h"
#include "math.h"
#include "sg.h"
#include "sg_mphase.h"
#include "flow.h"
#include "mem.h"
#include "metric.h"
#define R 8.314
#define pi 3.14159
#define T_SAT 373.15
#define h_lv 2455.1345e3
#define M 18.0
#define Beta 0.6666
#define a 23.224
#define b 3841.22
#define cc  -45.0

DEFINE_ADJUST(store_gradient, domain)
{
Thread *t;
Thread **pt;
cell_t c;
int phase_domain_index = 0.;
Domain *pDomain = DOMAIN_SUB_DOMAIN(domain,phase_domain_index);
Alloc_Storage_Vars(pDomain,SV_VOF_RG,SV_VOF_G,SV_NULL);
Scalar_Reconstruction(pDomain, SV_VOF,-1,SV_VOF_RG,NULL);
Scalar_Derivatives(pDomain,SV_VOF,-1,SV_VOF_G,SV_VOF_RG,Vof_Deriv_Accumulate);
mp_thread_loop_c (t,domain,pt)
if (FLUID_THREAD_P(t))
{
Thread *ppt = pt[phase_domain_index];

begin_c_loop (c,t)
{
C_UDMI(c,t,0) = C_VOF_G(c,ppt)[0];
C_UDMI(c,t,1) = C_VOF_G(c,ppt)[1];
C_UDMI(c,t,2) = NV_MAG(C_VOF_G(c,ppt));

}
end_c_loop (c,t)
}
Free_Storage_Vars(pDomain,SV_VOF_RG,SV_VOF_G,SV_NULL);
}
DEFINE_SOURCE(vap_src, cell, pri_th, dS, eqn)
{
Thread *mix_th, *sec_th;
real m_dot_v,T_v,T_l,P_v,P_sat;
mix_th = THREAD_SUPER_THREAD(pri_th);
sec_th = THREAD_SUB_THREAD(mix_th,1);
T_v=C_T(cell,pri_th);
P_v=C_P(cell,pri_th);
T_l=C_T(cell,sec_th);
P_sat=exp(a-(b/(T_l+cc)));
if(C_T(cell,mix_th)>=T_SAT)
{
m_dot_v = 0.0;

dS[eqn] = 0.0;
}

if(C_T(cell,mix_th)<=T_SAT)
{
m_dot_v =C_VOF(cell,sec_th)*(2.0*Beta/(2.0-Beta))*C_UDMI(cell,mix_th,2)*sqrt(M/(2*pi*R))*((P_v/sqrt(T_v))-(P_sat/sqrt(T_l)));
C_UDMI(cell,mix_th,3) =m_dot_v;
dS[eqn] = 0.0;
}
return m_dot_v ;
}
DEFINE_SOURCE(liq_src,cell,sec_th,dS,eqn)
{
Thread * mix_th, *pri_th;
real m_dot_l,T_v,T_l,P_v,P_sat;
mix_th=THREAD_SUPER_THREAD(sec_th);
pri_th=THREAD_SUB_THREAD(mix_th,0);
T_v = C_T(cell,pri_th);
P_v = C_P(cell,pri_th);
T_l = C_T(cell,sec_th);
P_sat=exp(a-(b/(T_l+cc)));
if(C_T(cell,mix_th)>=T_SAT)
{
m_dot_l = 0.0;
dS[eqn] = 0.0;
}
if(C_T(cell,mix_th)<=T_SAT)
{
m_dot_l = C_VOF(cell,sec_th)*(2.0*Beta/(2.0-Beta))*C_UDMI(cell,mix_th,2)*sqrt(M/(2*pi*R))*((P_v/sqrt(T_v))-(P_sat/sqrt(T_l)));
C_UDMI(cell,mix_th,4) =m_dot_l;
dS[eqn] =(2.0*Beta/(2.0-Beta))*C_UDMI(cell,mix_th,2)*sqrt(M/(2*pi*R))*((P_v/sqrt(T_v))-(P_sat/sqrt(T_l)));
}

return m_dot_l;
}
DEFINE_SOURCE(enrg_src,cell,mix_th,dS,eqn)
{
Thread *pri_th, *sec_th;
real m_dot,T_v,T_l,P_v,P_sat;
pri_th = THREAD_SUB_THREAD(mix_th,0);
sec_th = THREAD_SUB_THREAD(mix_th,1);
T_v=C_T(cell,pri_th);
P_v=C_P(cell,pri_th);
T_l=C_T(cell,sec_th);
P_sat=exp(a-(b/(T_l+cc)));
if(C_T(cell,mix_th)>=T_SAT)
{
m_dot =0.0;
dS[eqn] =0.0;
}

if(C_T(cell,mix_th)<=T_SAT)
{
m_dot=C_VOF(cell,sec_th)*(2.0*Beta/(2.0-Beta))*C_UDMI(cell,mix_th,2)*sqrt(M/(2*pi*R))*((P_v/sqrt(T_v))-(P_sat/sqrt(T_l)));
C_UDMI(cell,mix_th,5) =h_lv*m_dot;
dS[eqn]= 0.0;

}
return h_lv*m_dot; }

in this code I would change only 1 thing, which probably will not affect
was #define h_lv 2455.1345e3
to be #define h_lv 2.4551345e6

from your error I may assume you have problem with temperature convergence.
you may decrease time step
check your mesh
check boundary conditions

[pakk]

I tried to refactor your code, because you repeated your equations a few times:

#include "udf.h"
#include "math.h"
#include "sg.h"
#include "sg_mphase.h"
#include "flow.h"
#include "mem.h"
#include "metric.h"
#define R 8.314
#define pi 3.14159
#define T_SAT 373.15
#define h_lv 2455.1345e3
#define M 18.0
#define Beta 0.6666
#define a 23.224
#define b 3841.22
#define cc  -45.0

DEFINE_ADJUST(store_gradient, domain)
{
 Thread *t;
 Thread **pt;
 cell_t c;
 int phase_domain_index = 0.;
 Domain *pDomain = DOMAIN_SUB_DOMAIN(domain,phase_domain_index);
 Alloc_Storage_Vars(pDomain,SV_VOF_RG,SV_VOF_G,SV_NULL);
 Scalar_Reconstruction(pDomain, SV_VOF,-1,SV_VOF_RG,NULL);
 Scalar_Derivatives(pDomain,SV_VOF,-1,SV_VOF_G,SV_VOF_RG,Vof_Deriv_Accumulate);
 mp_thread_loop_c (t,domain,pt)
 if (FLUID_THREAD_P(t))
 {
  Thread *ppt = pt[phase_domain_index];
  begin_c_loop (c,t)
  {
   C_UDMI(c,t,0) = C_VOF_G(c,ppt)[0];
   C_UDMI(c,t,1) = C_VOF_G(c,ppt)[1];
   C_UDMI(c,t,2) = NV_MAG(C_VOF_G(c,ppt));
  }
 end_c_loop (c,t)
 }
Free_Storage_Vars(pDomain,SV_VOF_RG,SV_VOF_G,SV_NULL);
}

real ds(Thread* pri_th, Thread* sec_th, cell_t cell) {
 real T_v, T_l, P_v, P_sat, dS;
 if(C_T(cell,mix_th)>=T_SAT) {
  dS = 0;
 } else {
  T_v=C_T(cell,pri_th);
  P_v=C_P(cell,pri_th);
  T_l=C_T(cell,sec_th);
  P_sat=exp(a-(b/(T_l+cc)));dS;
  dS = (2.0*Beta/(2.0-Beta))*C_UDMI(cell,mix_th,2)*sqrt(M/(2*pi*R))*((P_v/sqrt(T_v))-(P_sat/sqrt(T_l)));
 }
 return dS;
}

real mdot(Thread* pri_th, Thread* sec_th, cell_t cell) {
 real T_v, T_l, P_v, P_sat, m_dot;
 if(C_T(cell,mix_th)>=T_SAT) {
  m_dot = 0;
 } else {
  T_v=C_T(cell,pri_th);
  P_v=C_P(cell,pri_th);
  T_l=C_T(cell,sec_th);
  P_sat=exp(a-(b/(T_l+cc)));
  m_dot = C_VOF(cell,sec_th)*(2.0*Beta/(2.0-Beta))*C_UDMI(cell,mix_th,2)*sqrt(M/(2*pi*R))*((P_v/sqrt(T_v))-(P_sat/sqrt(T_l)));
 }
 return m_dot;
}

DEFINE_SOURCE(vap_src, cell, pri_th, dS, eqn)
{
 Thread* mix_th = THREAD_SUPER_THREAD(pri_th);
 Thread* sec_th = THREAD_SUB_THREAD(mix_th,1);
 real m_dot_v=mdot(pri_th, sec_th, cell);
 C_UDMI(cell,mix_th,3) = m_dot_v;
 dS[eqn] = 0.0;
 return m_dot_v ;
}

DEFINE_SOURCE(liq_src,cell,sec_th,dS,eqn)
{
 Thread* mix_th=THREAD_SUPER_THREAD(sec_th);
 Thread* pri_th=THREAD_SUB_THREAD(mix_th,0);
 real m_dot_l=mdot(pri_th, sec_th, cell);
 C_UDMI(cell,mix_th,4) = m_dot_l;
 dS[eqn] = ds(pri_th, sec_th, cell);
 return m_dot_l;
}

DEFINE_SOURCE(enrg_src,cell,mix_th,dS,eqn)
{
 Thread* pri_th = THREAD_SUB_THREAD(mix_th,0);
 Thread* sec_th = THREAD_SUB_THREAD(mix_th,1);
 real m_dot=mdot(pri_th, sec_th, cell);
 C_UDMI(cell,mix_th,5) = h_lv*m_dot;
 dS[eqn]= 0.0;
 return h_lv*m_dot; 
}

This should be functionally the same (but I did not test it). What is now more apparent is that your two mass sources have the same sign, where I would expect them to differ in sign (whatever leaves the liquid zone should enter the vapor zone and vice versa).

I also noticed that only one of your sources specifies dS (as non-zero); the code above is (relatively) easy to change if you want to improve that.

[m.r.soufivand]

Thanks so much for your response and help. I try to use your code but gives the error is that mix_th is not defined.

[pakk]

I see, the definitions for ds and mdot should also have this:

Thread* mix_th = THREAD_SUPER_THREAD(pri_th);