Hi Rob,
Thank you for the response. I am not modelling the coal as a simple solid cell zone only. I am modelling the coal stockpile as a porous fluid zone with non-equilibrium thermal model enabled.
My intention is to reproduce a low-temperature coal self-heating/spontaneous oxidation model similar to Zhang et al. The physical idea is that air flows through the porous coal pile, oxygen is consumed by coal oxidation, CO₂ and CO are generated, and the reaction heat is deposited into the coal solid phase.
My current Fluent setup is as follows.
1. Cell zones
I have two main domains:
- Farfield air domain
- Cell zone name:
farfield - Type: fluid
- Zone ID: 6
- This represents the surrounding wind/air region.
- Coal stockpile domain
- Cell zone name:
stockpile - Type: fluid
- Zone ID: 7
- This is assigned as a porous fluid zone.
- The material assigned to this zone is the gas mixture used for species transport.
After enabling the Non-Equilibrium Thermal Model inside the porous zone, Fluent automatically created a corresponding solid zone:
- Phantom/solid coal zone
- Cell zone name:
stockpile:016 - Type: solid
- Zone ID: 16
- This is selected under the non-equilibrium thermal model as the solid zone coupled to the porous fluid zone.
So, the coal pile is not assigned as only a solid body. It is represented as a porous fluid region plus a coupled solid thermal zone through Fluent’s non-equilibrium porous model.
2. Models enabled
The following models are enabled:
- Energy equation: ON
- Species transport: ON
- Viscous model: RNG k-epsilon
- Multiphase: OFF
- Discrete phase: OFF
- Radiation: OFF
- Porous media: ON for the coal stockpile zone
- Non-equilibrium thermal model: ON for the stockpile porous zone
The reason for using the non-equilibrium thermal model is that I want separate gas and coal solid temperatures. The gas flows through the void spaces, while the coal matrix stores heat and receives heat generation from low-temperature oxidation.
3. Porous zone assignment
For the stockpile fluid zone, I enabled:
- Porous Zone
- Laminar Zone
- Source Terms
- Non-Equilibrium Thermal Model
The porous properties I assigned are:
The viscous resistance was calculated from the Carman–Kozeny permeability relation using:
This gives approximately:
For the non-equilibrium thermal settings, I used:
The interfacial area density is based on:
However, I am currently testing higher values such as 2400 1/m to strengthen gas-solid heat exchange and check sensitivity.
4. Material assignment
For the porous fluid zone, I assigned the gas mixture material.
The mixture contains:
The intended species order in the UDF is:
N₂ is treated as the background/bulk gas species.
For the coal solid zone, I assigned coal-like thermal properties:
5. Boundary conditions
The boundary conditions are:
Inlet
The inlet wind velocity is applied through a UDF profile using a 1/7 power-law wind profile:
where:
The inlet temperature is:
The inlet species mass fractions are:
Outlet
Walls
The bottom and top boundaries are assigned as walls. The bottom wall represents the ground, and the top wall is currently treated as a wall boundary following the paper-style simplified domain.
6. Initial condition and solution strategy
I am using a two-step solution method.
Step 1: Steady wind-flow solution
First, I solve the flow field without coal oxidation/heat generation. The goal is to obtain a stable wind field and pressure distribution around and inside the porous coal stockpile.
During this step, I keep the source terms either off or inactive.
Step 2: Transient self-heating simulation
Then I switch to transient simulation and activate the UDF source terms.
I patch/initialize the coal stockpile porous zone with:
I have tested both 293 K and a slightly elevated coal initial temperature around 303 K to provide a small thermal trigger for the Arrhenius oxidation reaction.
The transient time step is currently:
The reason is that low-temperature coal oxidation is a slow process and the paper also treats the heating process over days.
7. UDF source terms
I am using DEFINE_SOURCE UDFs for oxygen consumption, CO₂ generation, CO generation, total mass source, and heat generation.
The intended low-temperature oxidation reaction is:
The reaction rate is Arrhenius-type:
where:
The UDF source terms are intended as follows:
| UDF source | Hook location |
|---|
o2_sink | stockpile porous fluid zone → source term for O2 |
co2_source | stockpile porous fluid zone → source term for CO2 |
co_source | stockpile porous fluid zone → source term for CO |
mass_source | stockpile porous fluid zone → source term for total mass |
solid_heat | stockpile:016 solid zone → energy source |
wind_inlet | inlet velocity profile |
The heat source is applied to the solid coal zone, not directly to the gas phase, because the coal oxidation heat should be stored in the coal matrix first and then exchanged with the gas through the non-equilibrium thermal coupling.
8. Current issue
The UDF compiles successfully and I can hook the functions in Fluent. However, when I run the simulation, I do not see the expected temperature rise or hotspot formation. The result looks almost the same as if the UDF source terms are not active.
UDF:
#include "udf.h"
/* ===== Physical constants (Zhang Table 1) ===== */
#define A_PRE 1.8 /* Pre-exponential, 1/s */
#define EA 22000.0 /* Activation energy, J/mol */
#define R_UNIV 8.314 /* Gas constant, J/(mol K) */
#define DELTAH 3.0e8 /* Heat of reaction, J/kmol O2 */
#define EPS 0.2 /* Porosity */
/* Molar masses, kg/mol */
#define M_O2 0.032
#define M_CO2 0.044
#define M_CO 0.028
/* Wind profile */
#define V_REF 4.0
#define H_REF 10.0
/* Safety bounds */
#define T_MIN 250.0
#define T_MAX 1500.0
#define RHO_MIN 0.1
#define Y_FLOOR 0.0
#define Y_CEIL 1.0
#define Y_GROUND 1.0e-12
/* Debug: cell ID to print from (limits console spam) */
#define DEBUG_CELL 100
/* ============================================================
* Helper: volumetric reaction rate, kmol/m^3/s
* Reads SOLID temperature via THREAD_SUB_THREAD.
* Includes (1-eps) factor for solid volume fraction.
* ============================================================ */
real reaction_rate(cell_t c, Thread *t_fluid)
{
Thread *t_solid;
real T_coal, rho_gas, Y_O2, C_O2_kmol;
t_solid = THREAD_SUB_THREAD(t_fluid, 1);
if (t_solid == NULL)
return 0.0;
T_coal = C_T(c, t_solid);
if (T_coal < T_MIN) T_coal = T_MIN;
if (T_coal > T_MAX) T_coal = T_MAX;
rho_gas = C_R(c, t_fluid);
if (rho_gas < RHO_MIN) rho_gas = RHO_MIN;
Y_O2 = C_YI(c, t_fluid, 0);
if (Y_O2 < Y_FLOOR) Y_O2 = Y_FLOOR;
if (Y_O2 > Y_CEIL) Y_O2 = Y_CEIL;
C_O2_kmol = rho_gas * Y_O2 / M_O2 / 1000.0;
return (1.0 - EPS) * A_PRE * C_O2_kmol * exp(-EA / (R_UNIV * T_coal));
}
/* ============================================================
* Inlet wind profile: 1/7 power law
* ============================================================ */
DEFINE_PROFILE(wind_inlet, thread, position)
{
face_t f;
real x[ND_ND], y;
begin_f_loop(f, thread)
{
F_CENTROID(x, f, thread);
y = x[1];
if (y < Y_GROUND) y = Y_GROUND;
F_PROFILE(f, thread, position) = V_REF * pow(y / H_REF, 1.0/7.0);
}
end_f_loop(f, thread)
}
/* ============================================================
* O2 consumption (negative source)
* ============================================================ */
DEFINE_SOURCE(o2_sink, c, t, dS, eqn)
{
real r = reaction_rate(c, t);
real source_val = -r * M_O2 * 1000.0;
dS[eqn] = 0.0;
if (c == DEBUG_CELL)
Message("DEBUG o2_sink: cell=%d, r=%g, source=%g kg/m3/s\n",
c, r, source_val);
return source_val;
}
/* ============================================================
* CO2 production (positive source)
* ============================================================ */
DEFINE_SOURCE(co2_source, c, t, dS, eqn)
{
real r = reaction_rate(c, t);
real source_val = r * M_CO2 * 1000.0;
dS[eqn] = 0.0;
if (c == DEBUG_CELL)
Message("DEBUG co2_source: cell=%d, r=%g, source=%g kg/m3/s\n",
c, r, source_val);
return source_val;
}
/* ============================================================
* CO production (positive source, 0.1 stoichiometric ratio)
* ============================================================ */
DEFINE_SOURCE(co_source, c, t, dS, eqn)
{
real r = reaction_rate(c, t);
real source_val = 0.1 * r * M_CO * 1000.0;
dS[eqn] = 0.0;
if (c == DEBUG_CELL)
Message("DEBUG co_source: cell=%d, r=%g, source=%g kg/m3/s\n",
c, r, source_val);
return source_val;
}
/* ============================================================
* Heat generation on the SOLID coal phase
* ============================================================ */
DEFINE_SOURCE(solid_heat, c, t, dS, eqn)
{
Thread *t_fluid;
real r, source_val;
t_fluid = THREAD_SUPER_THREAD(t);
if (t_fluid == NULL)
return 0.0;
r = reaction_rate(c, t_fluid);
source_val = r * DELTAH;
dS[eqn] = 0.0;
if (c == DEBUG_CELL)
Message("DEBUG solid_heat: cell=%d, r=%g, source=%g W/m3\n",
c, r, source_val);
return source_val;
}
You are navigating away from the AIS Discovery experience