The basic equations of Darcy flow coupled to thermal transport are: $$$$\nabla \cdot \vec{u} = 0 \\ \vec{u} = -\frac{\mathbf{K}}{\mu} (\nabla p - \rho \vec{g}) \\ C\left( \frac{\partial T}{\partial t} + \epsilon \vec{u}\cdot\nabla T \right) - \nabla\cdot k \nabla T = 0$$$$ where $$$\vec{u}$$$ is the fluid velocity, $$$\epsilon$$$ is porosity, $$$\mathbf{K}$$$ is the permeability tensor, $$$\mu$$$ is fluid viscosity, $$$p$$$ is the pressure, $$$\vec{g}$$$ is the gravity vector, and $$$T$$$ is the temperature.

The parameters $$$\rho$$$, $$$C$$$, and $$$k$$$ are the porosity-dependent density, heat capacity, and thermal conductivity of the combined fluid/solid medium, defined by: $$$$\rho \equiv \epsilon \rho_f + (1-\epsilon) \rho_s \\ C \equiv \epsilon \rho_f {c_p}_f + (1-\epsilon) \rho_s {c_p}_s \\ k \equiv \epsilon k_f + (1-\epsilon) k_s$$$$ where $$$\epsilon$$$ is the porosity, $$$c_p$$$ is the specific heat, and the subscripts $$$f$$$ and $$$s$$$ refer to fluid and solid, respectively.

We shall henceforth assume that $$$\vec{g} = \vec{0}$$$. Taking the divergence of the second equation, and imposing the divergence-free condition leads to the following system of two equations in the unknowns $$$p$$$ and $$$T$$$: $$$$-\nabla \cdot \frac{\mathbf{K}}{\mu} \nabla p = 0 \\ C\left( \frac{\partial T}{\partial t} + \epsilon \vec{u}\cdot\nabla T \right) - \nabla \cdot k \nabla T = 0$$$$ where $$$\vec{u} = -\frac{\mathbf{K}}{\mu} \nabla p$$$. That is, the pressure satisfies Laplace's equation (with possibly non-constant coefficients) and the temperature satisfies the linear convection-diffusion equation.

We next compute some representative values for the mean fluid velocity induced by a given pressure head for the flow of water through a packed spherical medium, using the experimental apparatus described by Pamuk and Ozdemir (2012). Based on Fig. 2, we estimate the following values for the non-dimensionalized pressure drop $$$$P' \equiv \frac{\Delta P}{L} \frac{D^2}{\mu u_m}$$$$ and $$$$Re_D \equiv \frac{\rho u_m D}{\mu}$$$$

The Reynolds number based on the test chamber diameter:

Name | $$$P'$$$ | $$$Re_D$$$ |
---|---|---|

First Medium | $$$4 \times 10^6$$$ | 700 |

Second Medium | $$$0.5 \times 10^6$$$ | 950 |

Given these values and following geometric specifications and fluid properties of water at 30C:

Property | Value |
---|---|

Viscosity, $$$\mu$$$ | $$$7.98 \times 10^{-4}$$$ Pa-s |

Density, $$$\rho$$$ | $$$995.7$$$ kg/m$$$^3$$$ |

Test chamber length, $$$L$$$ | $$$304$$$ mm |

Test chamber diameter, $$$D$$$ | $$$51.4$$$ mm |

We obtain the following dimensional pressure drop and mean velocity values:

Name | $$$\Delta P$$$ | $$$u_m$$$ |
---|---|---|

First Medium | $$$4008.8$$$ Pa | $$$1.09 \times 10^{-2}$$$ m/s |

Second Medium | $$$680.07$$$ Pa | $$$1.48 \times 10^{-2}$$$ m/s |

The continuous Galerkin finite element method is known to exhibit spurious numerical oscillations unless the cell Peclet number, $$$Pe$$$, satisfies: $$$$Pe \equiv \frac{|\vec{u}| h}{2 \kappa} < 1$$$$ where $$$h$$$ is a representative finite element grid cell size, and $$$\kappa\equiv\frac{k}{\rho c}$$$ is the thermal diffusivity. All of the parameters in the definition of $$$Pe$$$ depend on the fluid/solid medium except the velocity magnitude $$$|\vec{u}|$$$, and the grid size $$$h$$$. This formula can therefore be used to determine a maximum-allowable cell size $$$h$$$ depending on the flow velocity magnitude.

The following thermal properties for water and steel were utilized.

Water | Value |
---|---|

Thermal conductivity, $$$k$$$ | $$$0.6$$$ W/m-K |

Specific heat capacity, $$$c_p$$$ | $$$4181.3$$$ J/kg-K |

Density, $$$\rho$$$ | $$$995.7$$$ kg/m$$$^3$$$ |

Steel | Value |
---|---|

Thermal conductivity, $$$k$$$ | $$$18$$$ W/m-K |

Specific heat capacity, $$$c_p$$$ | $$$466$$$ J/kg-K |

Density, $$$\rho$$$ | $$$8000$$$ kg/m$$$^3$$$ |

Using these properties and a porosity of $$$\epsilon = .25952$$$ yields the combined fluid/solid properties:

Combined | Value |
---|---|

Thermal conductivity, $$$k$$$ | $$$13.48$$$ W/m-K |

Specific heat capacity, $$$c_p$$$ | $$$1430.19$$$ J/kg-K |

Density, $$$\rho$$$ | $$$6182.24$$$ kg/m$$$^3$$$ |

and a combined thermal diffusivity of $$$$\kappa \approx 1.525 \times 10^{-6} \; \mathrm{m}^2/\mathrm{s}$$$$

Assuming a velocity of approximately 1 cm/s, the Peclet condition gives us the restriction $$$$h < \frac{2 \kappa}{|\vec{u}|} \approx 3 \times 10^{-4} \; \mathrm{m}$$$$

The cell Peclet number is only a necessary (but not sufficient) condition for numerical oscillations to be present in the solution, and is valid for the steady convection-diffusion equation. Numerical timestepping schemes may add artificial diffusion of $$$O(\Delta t)$$$, inhibiting numerical oscillations, but producing overly-diffusive solutions.