CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   Fluent UDF and Scheme Programming (https://www.cfd-online.com/Forums/fluent-udf/)
-   -   UDS of concentration solved in multiphases system (https://www.cfd-online.com/Forums/fluent-udf/182816-uds-concentration-solved-multiphases-system.html)

Tleja January 19, 2017 11:38

UDS of concentration solved in multiphases system
 
Hi Everyone!

I have done using UDS for multiphases system. I'm confused about Fluent solving. i 'm studying gas-solid reaction in fluidized bed.

I defined User-defined scalar as molar concentration (kmol/m3) in solid phase (secondary phase).
I follow equation 8–5 as shown in this link https://www.sharcnet.ca/Software/Ans..._uds_mp_theory
I neglected the diffusion term because of low impact to scalar.

Source term of UDS is a function of solid fraction. It means that where solid fraction appears high is where high value of source term is.

I could compile and calculate this code well, but i found weird result which UDS solved appears high at where low solid fraction is.

I also checked source term of the UDS. I found that source term appear high at where solid fraction is. It is reasonable source term. Basically, if i do right, UDS should increase with source term.

Could anyone help me to fix this problem?

I have faced this weird result for long time.

I already made consistency unit for all term (unsteady, convection and source term) in the unit of kg/m3s kmol/m3 =kg/m3s(scalar unit) as shown in coding below:

Code:

#include "udf.h"

DEFINE_SOURCE(UDS_con,c,t,dS,eqn)
{
        Thread *tm=THREAD_SUPER_THREAD(t);/*obtain mixture thread */
        Thread **pt = THREAD_SUB_THREADS(tm); /* pointer to sub_threads */
        Thread *tp = pt[0];                        /* Thread for primary phase */
        Thread *ts = pt[1];                        /* Thread for secondary phase*/
        real K,source;
       
        K=702000*exp(-40230/(8.314*C_T(c,ts)));      /* unit 1/s*/
       
       
        source=C_R(c,ts)*C_VOF(c,ts)*K*C_UDSI(c,ts,0);
        /*kg/m3s kmol/m3 =kg/m3s(scalar unit)*/

        dS[eqn]=C_R(c,ts)*C_VOF(c,ts)*K;
       
        return source;
       
}

DEFINE_UDS_UNSTEADY(unsteady,c,t,i,apu,su)
 {
    Thread *tm = THREAD_SUPER_THREAD(t); /*obtain mixture thread */
    Thread **pt = THREAD_SUB_THREADS(tm);  /* pointer to sub_threads */
    Thread *tp = pt[0];                        /* Thread for primary phase */
    Thread *ts = pt[1];    /* secondary phase thread */
    real physical_dt, vol, rho_m1, rho, vofs_m1, vofs, phi_old;
    physical_dt = RP_Get_Real("physical-time-step");
    vol = C_VOLUME(c,t);
    rho_m1 = C_R_M1(c,ts);
    rho = C_R(c,ts);
    vofs_m1 = C_VOF_M1(c,ts);
    vofs = C_VOF(c,ts);
    *apu = -rho*vol*vofs/physical_dt;/*implicit part*/
    phi_old = C_STORAGE_R(c,ts,SV_UDSI_M1(0));
    *su = rho_m1*vol*vofs_m1*phi_old/physical_dt;/*explicit part*/
 }

DEFINE_UDS_FLUX(my_uds_flux,f,t,i)
 {
        Thread *tm = THREAD_SUPER_THREAD(t); /*obtain mixture thread */
    Thread **pt = THREAD_SUB_THREADS(tm);  /* pointer to sub_threads */
    Thread *tp = pt[0];                        /* Thread for primary phase */
    Thread *ts = pt[1];    /* secondary phase thread */

    cell_t c0, c1 = -1;
    Thread *t0, *t1 = NULL;
    real NV_VEC(psi_vec), NV_VEC(A), flux = 0.0;
    c0 = F_C0(f,ts);
    t0 = F_C0_THREAD(f,ts);
    F_AREA(A, f, t);
    /* If face lies at domain boundary, use face values; */
    /* If face lies IN the domain, use average of adjacent cells. */
  if (BOUNDARY_FACE_THREAD_P(t)) /*Most face values will be available*/
    {
      real dens,vofs;
      /* Depending on its BC, density and viscosity may not be set on face thread*/
      if (NNULLP(THREAD_STORAGE(t,SV_DENSITY)))
        dens = F_R(f,ts);  /* Set dens to face value if available */
      else
        dens = C_R(c0,t0); /* else, set dens to cell value */
          if (NNULLP(THREAD_STORAGE(t,SV_VOF)))
        vofs = F_VOF(f,ts);  /* Set vofs to face value if available */
      else
        vofs = C_VOF(c0,t0); /* else, set vofs to cell value */

      NV_DS(psi_vec, =, F_U(f,ts), F_V(f,ts), F_W(f,ts), *,dens*vofs);
      flux = NV_DOT(psi_vec, A); /* flux through Face */ 
  }
    else
  {
    c1 = F_C1(f,ts);  /* Get cell on other side of face */
    t1 = F_C1_THREAD(f,ts);
    NV_DS(psi_vec, =, C_U(c0,t0),C_V(c0,t0),C_W(c0,t0),*,C_R(c0,t0)*C_VOF(c0,t0));
    NV_DS(psi_vec, +=, C_U(c1,t1),C_V(c1,t1),C_W(c1,t1),*,C_R(c1,t1)*C_VOF(c1,t1));
    flux = NV_DOT(psi_vec, A)/4.0; /* Average flux through face */
  }
       
  /* ANSYS Fluent will multiply the returned value by phi_f (the scalar’s
  value at the face) to get the ‘‘complete’’ advective term. */
  return flux;
 }

If i can't explain clearly, i'm so sorry for my bad English writing.:)

Kindly Regard
Tleja


All times are GMT -4. The time now is 14:02.