For those interested in some code (see post 2) 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

For those interested in some code (see post 2) 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;

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.

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

Well, the code is actually at:

for Dos/Windows: www.perkins-software.com.au/Fluent/udf_scalar_dos.c for Unix: www.perkins-software.com.au/Fluent/udf_scalar_unix.c

Have fun!

Greg

Greg,

a very interesting and helpful post. I just wanted to ask where did you find it particularly difficult to implement variable porosity in F5? I ask this because I've done this before (my PhD thesis, due to be handed in in a month, is based around this concept) and it wasn't as difficult. Then again, I've never used F4.5 since I've moved to Fluent straight from FIDAP, so maybe there is a big difference in implementation.

Anyway, this is just my curiosity since maybe I could have developed one third of my thesis in a third of the time! Go figure

Cheerio! AL

-- Alfonso Ferrandez | ferrandez@mech-eng.leeds.ac.uk Computational BioFluid Mechanics Research Group The University of Leeds, Leeds, UK URL: http://www.mech-eng.leeds.ac.uk/ferrandez/

Alfonso,

that's very intresting! How did you handle the change in porosity on the convective term in the various momentum and energy equations??????

When I looked into it, I worked out these terms need to be scaled by the porosity, which is a local cell property. Ie., similar to the equations for multiphase flow in Fluent v4.4/4.5 chapter 19.

One option was to correct for this with a source term, but this seemed very difficult to implement.

What did you do???? How???? Do you have the full set of equations you solve for in Fluent v5?? what type of UDFs did you write - different from mine???

Very interested to find out more,

Regards

Greg Perkins

Hi Greg!

That was a quick reply if I ever saw one

I did in fact write my own source term equations, since the porous media set in Fluent wasn't complex enough. This required a lot of literature searching on porous transport until eventually I came up with the proper source term functions.

Unfortunately I can't send you the details as of yet, since my thesis has been set under a 2/3 year NDA (industry sponsors, go figure!)

I had the same aversion to source term development as you, but after weeks of frustration, I made the bloody thing work! I suppose if I'd known about F4.5 I wouldn't have bothered

Cheers! AL

Greg,

I noticed your UDF deals with commbustion/gasification of carbon (char?) in a packed bed. Are you aware of the work of Prof. William Hallett at Univ. of Ottawa in Canada?

http://by.genie.uottawa.ca/profs/hallett/packedbd.html

Regards,

Henrik

Thanks Henrik,

yes I'm solving for underground coal gasification - in which part of the process can be represented like a packed bed surface gasifier....

Regards

Greg

Thanks Greg, very helpful code, do you mind if I ask you some questions regarding the scalar implementation? Regards Arturo