# Governing Equations

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. 

[](---)
# Values for Spherical Media

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)](http://www.sciencedirect.com/science/article/pii/S0894177711002640).  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 Cell Peclet Restriction

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.