[ANSYS, Inc. Logo] return to home search
next up previous contents index

8.2.5 User-Defined Scalars

This section contains examples of UDFs that can be used to customize user-defined scalar (UDS) transport equations. See Section  2.7 in the UDF Manual for information on how you can define UDFs in ANSYS FLUENT. See this section in the separate Theory Guide for UDS equation theory and details on how to set up scalar equations.



Postprocessing Using User-Defined Scalars

Below is an example of a compiled UDF that computes the gradient of temperature to the fourth power, and stores its magnitude in a user-defined scalar. The computed temperature gradient can, for example, be subsequently used to plot contours. Although the practical application of this UDF is questionable, its purpose here is to show the methodology of computing gradients of arbitrary quantities that can be used for postprocessing. 

/***********************************************************************/
/* UDF for computing the magnitude of the gradient of T^4              */
/***********************************************************************/

#include "udf.h"

/* Define which user-defined scalars to use. */
enum
{
  T4,
  MAG_GRAD_T4,
  N_REQUIRED_UDS
};

DEFINE_ADJUST(adjust_fcn, domain)
{
  Thread *t;
  cell_t c;
  face_t f;


  /* Make sure there are enough user-defined scalars. */
  if (n_uds < N_REQUIRED_UDS)
    Internal_Error("not enough user-defined scalars allocated");

  /* Fill first UDS with temperature raised to fourth power. */  
  thread_loop_c (t,domain)
    {
      if (NULL != THREAD_STORAGE(t,SV_UDS_I(T4)))
        {
          begin_c_loop (c,t)
            {
              real T = C_T(c,t);
              C_UDSI(c,t,T4) = pow(T,4.);
            }
          end_c_loop (c,t)
        }
    }


  thread_loop_f (t,domain)
    {
      if (NULL != THREAD_STORAGE(t,SV_UDS_I(T4)))
        {
          begin_f_loop (f,t)
            {
              real T = 0.;
              if (NULL != THREAD_STORAGE(t,SV_T))
                T = F_T(f,t);
              else if (NULL != THREAD_STORAGE(t->t0,SV_T))
                T = C_T(F_C0(f,t),t->t0);
              F_UDSI(f,t,T4) = pow(T,4.);
            }
          end_f_loop (f,t)
        }
    }

  /* Fill second UDS with magnitude of gradient. */
  thread_loop_c (t,domain)
    {
      if (NULL != THREAD_STORAGE(t,SV_UDS_I(T4)) &&
          NULL != T_STORAGE_R_NV(t,SV_UDSI_G(T4)))
       {
          begin_c_loop (c,t)
            {
              C_UDSI(c,t,MAG_GRAD_T4) = NV_MAG(C_UDSI_G(c,t,T4));
            }
          end_c_loop (c,t)
        }
    }

  thread_loop_f (t,domain)
    {
      if (NULL != THREAD_STORAGE(t,SV_UDS_I(T4)) &&
          NULL != T_STORAGE_R_NV(t->t0,SV_UDSI_G(T4)))
       {
          begin_f_loop (f,t)
           {
            F_UDSI(f,t,MAG_GRAD_T4)=C_UDSI(F_C0(f,t),t->t0,MAG_GRAD_T4);
           }
          end_f_loop (f,t)
        }
    }
}

The conditional statement if (NULL != THREAD_STORAGE(t,SV_UDS_I(T4))) is used to check if the storage for the user-defined scalar with index T4 has been allocated, while NULL != T_STORAGE_R_NV(t,SV_UDSI_G(T4)) checks whether the storage of the gradient of the user-defined scalar with index T4 has been allocated.

In addition to compiling this UDF, as described in Chapter  5, you will need to enable the solution of a user-defined scalar transport equation in ANSYS FLUENT.

Define $\rightarrow$ User-Defined $\rightarrow$ Scalars...

See this section in the separate Theory Guide for UDS equation theory and details on how to setup scalar equations.



Implementing ANSYS FLUENT's P-1 Radiation Model Using User-Defined Scalars

This section provides an example that demonstrates how the P1 radiation model can be implemented as a UDF, utilizing a user-defined scalar transport equation. In the P1 model, the variation of the incident radiation , $G$, in the domain can be described by an equation that consists of a diffusion and source term.

The transport equation for incident radiation, $G$, is given by Equation  8.2-5. The diffusion coefficient, $\Gamma$, is given by Equation  8.2-6 and the source term is given by Equation  8.2-7. See this section in the separate Theory Guide for more details.


\nabla \cdot \left( {\Gamma} \nabla G \right) + S^G = 0 (8.2-5)


{\Gamma} = \frac{1}{3 a + \left( 3 - C \right) \sigma_s } (8.2-6)


S^G = a \left( 4\sigma T^4 - G \right) (8.2-7)

The boundary condition for $G$ at the walls is equal to the negative of the radiative wall heat flux, $q_{r,w}$ (Equation  8.2-8), where $\vec{n}$ is the outward normal vector (see this section in the separate Theory Guide for more details). The radiative wall heat flux can be given by Equation  8.2-9.


q_r \cdot {\vec n} = - {\Gamma} \nabla G \cdot {\vec n} (8.2-8)


q_{r,w} = - \frac{\epsilon_w}{2 \left(2 - \epsilon_w \right)} \left( 4 \sigma T_w^4 - G_w \right) (8.2-9)

This form of the boundary condition is unfortunately specified in terms of the incident radiation at the wall, $G_w$. This mixed boundary condition can be avoided by solving first for $G_w$ using Equations  8.2-8 and 8.2-9, resulting in Equation  8.2-10. Then, this expression for $G_w$ is substituted back into Equation  8.2-9 to give the radiative wall heat flux $q_{r,w}$ as Equation  8.2-11.


G_w = \frac{4 \sigma T_w^4 E_w + \frac{\alpha_0 {\Gamma}_0}{... ..._0 - \beta_0(G) \right]}{E_w + \frac{\alpha_0 {\Gamma}_0}{A}} (8.2-10)


q_{r} = - \; \frac{\alpha_0 {\Gamma}_0 E_w}{A \left( E_w + \... ...A} \right)} \left[4 \pi I_b(T_{iw}) -G_0 + \beta_0(G) \right] (8.2-11)

The additional $\beta_0$ and $G_0$ terms that appear in Equations  8.2-10 and 8.2-11 are a result of the evaluation of the gradient of incident radiation in Equation  8.2-8.

In ANSYS FLUENT, the component of a gradient of a scalar directed normal to a cell boundary (face), ${\nabla}G{\cdot}$ n, is estimated as the sum of primary and secondary components. The primary component represents the gradient in the direction defined by the cell centroids, and the secondary component is in the direction along the face separating the two cells. From this information, the face normal component can be determined. The secondary component of the gradient can be found using the ANSYS FLUENT macro BOUNDARY_SECONDARY_GRADIENT_SOURCE (which is described in Section  3.2.5). The use of this macro first requires that cell geometry information be defined, which can be readily obtained by the use of a second macro, BOUNDARY_FACE_GEOMETRY (see Section  3.2.5). You will see these macros called in the UDF that defines the wall boundary condition for $G$.

To complete the implementation of the P1 model, the radiation energy equation must be coupled with the thermal energy equation. This is accomplished by modifying the source term and wall boundary condition of the energy equation. Consider first how the energy equation source term must be modified. The gradient of the incident radiation is proportional to the radiative heat flux. A local increase (or decrease) in the radiative heat flux is attributable to a local decrease (or increase) in thermal energy via the absorption and emission mechanisms. The gradient of the radiative heat flux is therefore a (negative) source of thermal energy. The source term for the incident radiation Equation  8.2-7 is equal to the gradient of the radiative heat flux and hence its negative specifies the source term needed to modify the energy equation (see this section in the separate Theory Guide for more details).

Now consider how the energy boundary condition at the wall must be modified. Locally, the only mode of energy transfer from the wall to the fluid that is accounted for by default is conduction. With the inclusion of radiation effects, radiative heat transfer to and from the wall must also be accounted for. (This is done automatically if you use ANSYS FLUENT's built-in P1 model.) The DEFINE_HEAT_FLUX macro allows the wall boundary condition to be modified to accommodate this second mode of heat transfer by specifying the coefficients of the $qir$ equation discussed in Section  2.3.9. The net radiative heat flux to the wall has already been given as Equation  8.2-9. Comparing this equation with that for $qir$ in Section  2.3.9 will result in the proper coefficients for $cir[]$.

In this example, the implementation of the P1 model can be accomplished through six separate UDFs. They are all included in a single source file, which can be executed as a compiled UDF. The single user-defined scalar transport equation for incident radiation, $G$, uses a DEFINE_DIFFUSIVITY UDF to define $\Gamma$ of Equation  8.2-6, and a UDF to define the source term of Equation  8.2-7. The boundary condition for $G$ at the walls is handled by assigning, in DEFINE_PROFILE, the negative of Equation  8.2-11 as the specified flux. A DEFINE_ADJUST UDF is used to instruct ANSYS FLUENT to check that the proper number of user-defined scalars has been defined (in the solver). Lastly, the energy equation must be assigned a source term equal to the negative of that used in the incident radiation equation and the DEFINE_HEAT_FLUX UDF is used to alter the boundary conditions at the walls for the energy equation.

In the solver, at least one user-defined scalar (UDS) equation must be enabled. The scalar diffusivity is assigned in the Create/Edit Materials dialog box for the scalar equation. The scalar source and energy source terms are assigned in the boundary condition dialog box for the fluid zones. The boundary condition for the scalar equation at the walls is assigned in the boundary condition dialog box for the wall zones. The DEFINE_ADJUST and DEFINE_HEAT_FLUX functions are assigned in the User-Defined Function Hooks dialog box.

Note that the residual monitor for the UDS equation should be reduced from $1e-3$ to $1e-6$ before running the solution. If the solution diverges, then it may be due to the large source terms. In this case, the under-relaxation factor should be reduced to $0.99$ and the solution re-run. 

/**************************************************************/
/* Implementation of the P1 model using user-defined scalars  */
/**************************************************************/

#include "udf.h"
#include "sg.h"

/* Define which user-defined scalars to use. */
enum
{
  P1,
  N_REQUIRED_UDS
};


static real abs_coeff = 0.2;    /* absorption coefficient */
static real scat_coeff = 0.0;   /* scattering coefficient */
static real las_coeff = 0.0;    /* linear-anisotropic     */
                                /* scattering coefficient */
static real epsilon_w = 1.0;    /* wall emissivity */


DEFINE_ADJUST(p1_adjust, domain)
{
  /* Make sure there are enough user defined-scalars. */
  if (n_uds < N_REQUIRED_UDS)
    Internal_Error("not enough user-defined scalars allocated");
}


DEFINE_SOURCE(energy_source, c, t, dS, eqn)
{
  dS[eqn] = -16.*abs_coeff*SIGMA_SBC*pow(C_T(c,t),3.);
  return -abs_coeff*(4.*SIGMA_SBC*pow(C_T(c,t),4.) - C_UDSI(c,t,P1));
}


DEFINE_SOURCE(p1_source, c, t, dS, eqn)
{
  dS[eqn] = 16.*abs_coeff*SIGMA_SBC*pow(C_T(c,t),3.);
  return abs_coeff*(4.*SIGMA_SBC*pow(C_T(c,t),4.) - C_UDSI(c,t,P1));
}


DEFINE_DIFFUSIVITY(p1_diffusivity, c, t, i)
{
  return 1./(3.*abs_coeff + (3. - las_coeff)*scat_coeff);
}


DEFINE_PROFILE(p1_bc, thread, position)

{
  face_t f;
  real A[ND_ND],At;
  real dG[ND_ND],dr0[ND_ND],es[ND_ND],ds,A_by_es;
  real aterm,alpha0,beta0,gamma0,Gsource,Ibw;
  real Ew = epsilon_w/(2.*(2. - epsilon_w));
  Thread *t0=thread->t0;

  /* Do nothing if areas aren't computed yet or not next to fluid. */
  if (!Data_Valid_P() || !FLUID_THREAD_P(t0)) return;

  begin_f_loop (f,thread)
    {
      cell_t c0 = F_C0(f,thread);

      BOUNDARY_FACE_GEOMETRY(f,thread,A,ds,es,A_by_es,dr0);
      At = NV_MAG(A);

      if (NULLP(T_STORAGE_R_NV(t0,SV_UDSI_G(P1))))
        Gsource = 0.;   /* if gradient not stored yet */
      else
        BOUNDARY_SECONDARY_GRADIENT_SOURCE(Gsource,SV_UDSI_G(P1),
                                           dG,es,A_by_es,1.);

      gamma0 = C_UDSI_DIFF(c0,t0,P1);
      alpha0 = A_by_es/ds;
      beta0  = Gsource/alpha0;
      aterm  = alpha0*gamma0/At;

      Ibw = SIGMA_SBC*pow(WALL_TEMP_OUTER(f,thread),4.)/M_PI;

      /* Specify the radiative heat flux. */
      F_PROFILE(f,thread,position) =
        aterm*Ew/(Ew + aterm)*(4.*M_PI*Ibw - C_UDSI(c0,t0,P1) + beta0);
    }
  end_f_loop (f,thread)
}


DEFINE_HEAT_FLUX(heat_flux, f, t, c0, t0, cid, cir)
{

     real Ew = epsilon_w/(2.*(2. - epsilon_w));
  
     cir[0] = Ew * F_UDSI(f,t,P1);
     cir[3] = 4.0 * Ew * SIGMA_SBC;

}

next up previous contents index Previous: 8.2.4 Reaction Rates
Up: 8.2 Detailed UDF Examples
Next: 8.2.6 User-Defined Real Gas
Release 12.0 © ANSYS, Inc. 2009-01-14