UDFs for Scalar Eqn - Fluid/Solid HT

Greg Perkins

For those interested in some code for scalar transport eqns under Fluent v5 here's some examples. . . I developed this for solving heat transfer between fluid and a packed bed solid (after a few weeks I moved to Fluent v4.5 which is better for this application).

Anyway this should allow you to solve for solid phase temperature as a scalar transport eqn as a special user porous media model assuming a packed bed of spheres. I've cut down the code so it may or may not compile and I don't provide any guarantees . . it may be more confusing than helpful!!!

But if you play around I think you can get it to work - once compiled - use makefiles - don't think it will work under interpreted mode (but try??) - then link the source terms to the appropriate Fluent eqn. That's about it.

Unfortunately, if it doesn't work/doesn't work for you I can't provide detailed responses to questions - sorry.

A bit on how it works: it uses two main UDS's - ZETA = local porosity and SOLID_TEMP for local solid temperature. ZETA just stores the value, SOLID_TEMP can be solved for. Another six UDS are required for other debug type results that I found useful. So in total declare space for 8 UDS, but only solve for SOLID_TEMP - turn off solution of all others. (with mods you can update ZETA due to reaction for example) I would use the UDMI function in 5.4.8 for storing values at cells - but that wasn't available until recently.

There are source terms for momentum - SRCE_Xmon etc. - these calculate the drag due to the solid spheres on the fluid phase using the Ergun eqn. Link these into the fluid domain sections in the Boundary Conditions panel.

There are source terms for energy - fluid and solid - this just calculates the interphase heat transfer with a HTC correlation. Link these as source terms for fluid energy as solved by Fluent and the SOLID_TEMP UDS. If you solve steady probs there's really no point for this.

Use the diffusivity routine to tell Fluent how to calculate the diffusivity for your UDS's. This only calculates for the SOLID_TEMP eqn where it uses alpha = k/rho*Cp = thermal diffusivity.

You can use the init routine to help initialise the domain but will need mods depending on what you try out. Then set some b/cs - note to set also your SOLID_TEMP UDS b/c on each boundary - otherwise your solid temp will be 0 Kelvin which is not realistic!!!

Greg Perkins

ps: if you want to use a variable porosity use Fluent v4.5 - its toooo hard to implement in Fluent v5.x

/************************************************** **********************/ /* Reaction Module Source Code */ /* */ /* By Greg Perkins */ /* CANCES, Australian Technology Park */ /* Ph: 02 9318 0004, Fax: 02 9319 2328 */ /* Email: perkinsg@cances.atp.com.au */ /* */ /* Started: 03-01-2000 */ /* */ /* Revisions: */ /* */ /* 14-02-2000

(i) Solid phase temperature and UDS

(ii) General species structure, material properties and reactions */ /************************************************** **********************/

/* ---- These routines implement a hetrogeneous reaction model for the

combustion and gasification reactions of carbon together with a

porous media model */

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

/* ------------------------ REACTION MODEL DEFINES --------------------- */

#define PI 3.141592654 #define KJCAL 4.184 #define UNIVERSAL_GAS_CONST 8.314 /* -- [J/mol-K] */

#define SOLID_DENSITY(T) 1200 /* -- Coal material density [kg/m3] */ #define SOLID_THERMAL_COND(T) 1.47/1.0*pow(T/273,0.5)

/* -- Hobbs et al. 1993, [W/m-K] */ #define SOLID_CP(T) 3000.0*UNIVERSAL_GAS_CONST/12.0*(exp(1200/T)/pow((exp(1200/T)-1)*T/1200,2.0))

/* -- Hobbs et al. 1993, [J/kg-K] */

#define alpha_s 0.5 /* -- fraction of energy released in solid phase */ #define alpha_f 0.5 /* -- fraction of energy released in fluid phase */

#define ZETA_CO 0.8 /* -- Initial coal vol. frac. */ #define ZETA_VO 0.2 /* -- Initial void vol. frac. */ #define R_CO 0.01 /* -- Initial coal particle radius [m] */

#define ZETA_MAX 0.95 /* -- Maximum porosity without checks */ #define DTHI_MAX -1.0e30

#define PERM_EXP_VOID 25.0 #define TEMP_KELVIN 0 #define TEMP_DEGC 1 #define TREF 273.0

/* -------------------- POROUS MEDIA MODEL DEFINES --------------------- */

/* ---- use a user defined scalar to track material properties of the domain */ enum { ZETA, SOLID_TEMP, SRCE_SOLID, SRCE_FLUID, DTEMP, DZETA, THI, RE, N_REQUIRED_UDS };

/* ------------------------------------------------------------------------

Material Properties

This routine returns the local properties of the material (C2,C2) for

a given cell in the domain.

1. CHECK IMPLEMENTATION

------------------------------------------------------------------------ */ void Material_Properties(cell_t c, Thread *t, real *C1, real *C2) {

real zeta, zeta3, temp1;

/* --- use this routine to return the permeability, alpha and the

co-efficient C2, for each location x,y,z in the domain. No

Modifications needed for 2D. For time dependent porous media

use the RP_Get_Real("flow time") function to find out t (secs) */

zeta = C_UDSI(c,t,ZETA);

zeta3 = pow(zeta,3.0);

temp1 = pow((1.0-((zeta-ZETA_VO)/ZETA_CO)),(1.0/3.0));

temp1 = 1.0; /* ---- Dp fixed */

/* ---- Calculate local permeability = 1/C1 */

*C1 = (150.0*pow((1.0-zeta),2.0))/(zeta3*4.0*pow(R_CO,2.0)*pow(temp1,2.0));

/* ---- Calculate local C2 co-efficient */

*C2 = 1.75*(1.0-zeta)/(R_CO*temp1*zeta3); }

/* ------------------------------------------------------------------------

X_Momentum_Source

This routine returns the source term for the X-momentum term for each

control volume in the domain. The local properties are obtained by

calling Material_Properties

------------------------------------------------------------------------ */ DEFINE_SOURCE(SRCE_Xmom,cell,thread,dS,eqn) {

real C1, C2, constant1, constant2, Ux, source;

/* --- determine local properties */

Material_Properties(cell, thread, &C1, &C2);

/* --- determine constants 1,2 */

constant1 = C_MU_L(cell,thread) * C1;

constant2 = 0.5 * C_R(cell,thread) * C2;

/* --- determine x-velocity */

Ux = C_U(cell,thread);

source = -(constant1*Ux + constant2 * fabs(Ux) * Ux);

dS[eqn] = -(constant1 + 2 * constant2 * fabs(Ux));

return source; }

/* ------------------------------------------------------------------------

Y_Momentum_Source

This routine returns the source term for the Y-momentum term for each

control volume in the domain. The local properties are obtained by

calling Material_Properties

------------------------------------------------------------------------ */ DEFINE_SOURCE(SRCE_Ymom, cell, thread, dS, eqn) {

real C1, C2, constant1, constant2, Uy, source;

/* --- determine local properties */

Material_Properties(cell, thread, &C1, &C2);

/* --- determine constants 1,2 */

constant1 = C_MU_L(cell,thread) * C1;

constant2 = 0.5 * C_R(cell,thread) * C2;

/* --- determine y-velocity */

Uy = C_V(cell,thread);

source = -(constant1*Uy + constant2 * fabs(Uy) * Uy);

dS[eqn] = -(constant1 + 2 * constant2 * fabs(Uy));

return source; }

/* ------------------------------------------------------------------------

Z_Momentum_Source

This routine returns the source term for the Z-momentum term for each

control volume in the domain. The local properties are obtained by

calling Material_Properties

------------------------------------------------------------------------ */ DEFINE_SOURCE(SRCE_Zmom, cell, thread, dS, eqn) {

real C1, C2, constant1, constant2, Uz, source;

/* --- determine local properties */

Material_Properties(cell, thread, &C1, &C2);

/* --- determine constants 1,2 */

constant1 = C_MU_L(cell,thread) * C1;

constant2 = 0.5 * C_R(cell,thread) * C2;

/* --- determine z-velocity */

Uz = C_W(cell,thread);

source = -(constant1*Uz + constant2 * fabs(Uz) * Uz);

dS[eqn] = -(constant1 + 2 * constant2 * fabs(Uz));

return source; }

/* ------------------------------------------------------------------------

Initialization UDF

This routine initializes the solution in the computational domain.

------------------------------------------------------------------------ */ DEFINE_INIT(INIT_Solution, domain) {

cell_t c;

face_t f;

Thread *thread;

real x[ND_ND];

/* ---- Set initial porosity to be ZETA_VO for y>0.5m */

thread_loop_c (thread, domain)

{

begin_c_loop (c, thread)

{

C_CENTROID(x,c,thread);

if (x[1] > 0.0)

C_UDSI(c, thread, ZETA) = ZETA_VO;

else

C_UDSI(c, thread, ZETA) = 1.0;

}

end_c_loop (c, thread)

}

}

/* ------------------------------------------------------------------------

UDS Diffusivity [m2/s]

This routine returns the diffusivity for the user defined scalar

transport equations. The only non-zero setting is for the

solid phase temperature eqn (solid energy eqn). (RB-002, p101)

------------------------------------------------------------------------ */ DEFINE_DIFFUSIVITY(UDS_Diffusivity, cell, thread, i) {

real source,Ts;

if (i==SOLID_TEMP)

{

Ts = C_UDSI(cell,thread,SOLID_TEMP);

source = SOLID_THERMAL_COND(Ts)/( SOLID_DENSITY(Ts)*SOLID_CP(Ts) );

}

else source = 0.0;

return source; }

/* ------------------------------------------------------------------------

GAS-SOLID Heat Transfer Co-efficient [W/m2-K]

This routine calculates the gas-solid heat transfer co-efficient

source: Wakao & Kaguei, (1982), (RB-002, p103)

------------------------------------------------------------------------ */ real Gas_Solid_HTC(cell_t cell, Thread *thread) {

real Pr,Re,Us;

/* ---- calculate superficial velocity */

Us = sqrt( pow(C_U(cell,thread),2.0) + pow(C_V(cell,thread),2.0) +

pow(C_W(cell,thread),2.0) );

/* ---- calculate Prandtl number */

Pr = C_MU_EFF(cell,thread)*C_CP(cell,thread)/C_K_EFF(cell,thread);

/* ---- calculate Reynolds no. */

Re = 2*R_CO*Us*C_R(cell,thread)/C_MU_EFF(cell,thread);

return C_K_EFF(cell,thread)/(2*R_CO)*(2.0 + 1.1*pow(Pr,1.0/3.0)*pow(Re,0.6)); }

/* ------------------------------------------------------------------------

ENERGY SOURCE: Solid Phase [K/s]

This routine calculates the solid phase energy source term. This is

composed of two parts: the heat transfer between the solid and gas

and the energy released in the solid phase, (RB-002, p102)

------------------------------------------------------------------------ */ DEFINE_SOURCE(SRCE_Solid_Energy,cell,thread,dS,eqn ) {

real zeta, source, thi, dthi, Ts, HTC;

/* ---- Current porosity */

zeta = C_UDSI(cell, thread, ZETA);

/* ---- Specific Surface Area Term */

thi = (3*ZETA_CO/R_CO)*pow((1.0-((zeta-ZETA_VO)/ZETA_CO)),2.0/3.0);

if (zeta > ZETA_MAX)

dthi = DTHI_MAX;

else

dthi = -(2.0/R_CO)*pow((1.0-((zeta-ZETA_VO)/ZETA_CO)),-1.0/3.0);

/* --- Check that there is some solid in this cell*/

if (zeta>0.99) zeta=0.9999;

/* --- Solid Temperature & Heat Transfer Co-eff */

Ts = C_UDSI(cell,thread,SOLID_TEMP);

HTC = Gas_Solid_HTC(cell,thread);

source = -HTC*thi*(Ts-C_T(cell,thread));

source = source/( SOLID_DENSITY(Ts)*SOLID_CP(Ts)*(1-zeta) );

dS[eqn] = -HTC*thi/( SOLID_DENSITY(Ts)*SOLID_CP(Ts)*(1-zeta) );

C_UDSI(cell,thread,SRCE_SOLID) = source;

C_UDSI(cell,thread,RE) = HTC;

return source; }

/* ------------------------------------------------------------------------

ENERGY SOURCE: Fluid Phase [J/m3-s]

This routine calculates the fluid phase energy source term. This is

composed of two parts: the heat transfer between the solid and gas

and the energy released in the fluid phase, (RB-002, p102)

------------------------------------------------------------------------ */ DEFINE_SOURCE(SRCE_Fluid_Energy,cell,thread,dS,eqn ) {

real zeta, source, thi, dthi, Ts, HTC;

/* ---- Current porosity */

zeta = C_UDSI(cell, thread, ZETA);

/* ---- Specific Surface Area Term */

thi = (3*ZETA_CO/R_CO)*pow((1.0-((zeta-ZETA_VO)/ZETA_CO)),2.0/3.0);

if (zeta > ZETA_MAX)

dthi = DTHI_MAX;

else

dthi = -(2.0/R_CO)*pow((1.0-((zeta-ZETA_VO)/ZETA_CO)),-1.0/3.0);

/* --- Check that there is some fluid in this cell */

if (zeta<0.01) zeta=0.0001;

/* --- Solid Temperature & Heat Transfer Co-eff */

Ts = C_UDSI(cell,thread,SOLID_TEMP);

HTC = Gas_Solid_HTC(cell,thread);

source = HTC*thi*(Ts-C_T(cell,thread))/zeta;

dS[eqn] = -HTC*thi/zeta;

C_UDSI(cell,thread,SRCE_FLUID) = source;

C_UDSI(cell,thread,DTEMP) = (Ts - C_T(cell,thread) );

return source; }

/* -----------------------------------------------------------------------

REFERENCES

Felder, R. M., and Rousseau, R. W., (1986), "Elementary principles of

Chemical Processes",2nd Edition, John Wiley & Sons, New York.

Hobbs, M., L., Radulovic, P. T., and Smoot, L. D., (1993), "Combustion

and Gasification of Coals in Fixed Beds", Prog. Eenergy Combust. Sci.,

Vol. 19, pp.505-586.

Wakao, N., and Kaguei, S., (1982), "Heat and Mass Transfer in Packed

Beds", Gordon & Breach Science Publishers, New York.

----------------------------------------------------------------------- */