Darcy Flow - Solution by Finite Difference

Question

I'd like to model the fluid flow through a porous medium using finite differences. Since I am new to this numerical technique, I have a simple question.

I use the following set of equations to calculate the change in water volume or porosity over time.

$$q_x = K \frac{dP}{dx}$$

$$q_y = K \left(\frac{dP}{dy} + \varphi g\right)$$

$$\frac{d\phi}{dt} = \frac{q_x}{dx} + \frac{q_y}{dy} $$

where $K$, $P$, $g$, $q$, $\varphi$, $\phi$ are a constant, pressure, gravity, flow velocity, fluid density and porosity.

I know how to solve the groundwater flow equation for fluid flow along the x-axis:

$$\frac{d\phi_{i,j}}{dt} = K \frac{P_{i-1,j} - 2P_{i,j} + P_{i+1,j}}{dx^2 } + ??? $$

But how do I incorporate the gravity term for fluid flow along the y-axis? Would this be the correct solution? $$\frac{d\phi_{i,j}}{dt} = K \frac{P_{i-1,j} - 2P_{i,j} + P_{i+1,j}}{dx^2 } + K \frac{P_{i,j-1} - 2P_{i,j} + P_{i,j+1} + \varphi g}{dy^2} $$

Answer 1

No, that's not right.

First, I assume you mean partial derivatives everywhere you have written $d$. Second, your formula: $$ \frac{d \phi}{dt} = \frac{q_x}{dx} + \frac{q_y}{dy} $$ should probably be: $$ \frac{\partial \phi}{\partial t} = \frac{\partial q_x}{\partial x} + \frac{\partial q_y}{\partial y} $$ or $$ \frac{\partial \phi}{\partial t} = \nabla \cdot \boldsymbol{q} $$ Next, I'm assuming that where you write the gravity term with $\varphi$ you meant to write $\phi$, otherwise, you need to give a definition of $\varphi$.

Finally, given all of the above, you can see that your $y$-direction term should have a first derivative of $\phi$ in it: $$ \frac{\partial q_y}{\partial y} = \frac{\partial}{\partial y} \left(K\frac{\partial P}{\partial y} + \phi g\right) = \frac{\partial}{\partial y}\left(K \frac{\partial P}{\partial y}\right) + g\frac{\partial \phi}{\partial y} $$ Here I've assumed that $K$ was only supposed to multiply the first term in the definition of $q_y$, but if it's both, then you should have a $K$ in the coefficient of second term as well.

Based on this, you should have an approximation of the first-derivative of $\phi$ appearing in your formula. Given that this makes this problem a convection-diffusion problem, you may need to be very careful with how you implement that.

Answer 2

First up, all the things Bill said in his answer are spot on. You have to be fairly rigorous about the nomenclature here. The old-school expression you see for "Darcys Law" is something like the following:

$$ Q = \frac{-K A}{\mu} \frac{(P_b-P_a)}{L} $$

It is more often presented in the form given by

$$ \dot{Q} = - \frac{K}{\mu} \nabla P $$

which provides a simple stimulus response relationship analogous to Ohm's law in electronics and Fick's law of diffusion. In this presentation, the geometric terms (and gravity!) have been absorbed to produce a volumetric flow velocity $\dot{Q}$ in the direction of the gradient of pressure $\nabla P$, linked scalar viscosity and permeability of the media. Think of gravity as a force that creates a pressure gradient.

Darcy flow is inherently steady state and with a bit of hoop jumping you can get a pair of equations which are solvable to get satisfaction of the Darcy flow law and continuity. Omitting the math and symbol juggling, you rearrange things to get a pair of equations like this:

$$ \tilde{\nabla} \tilde{P}= \tilde{\nabla}^2 \tilde{v} $$

$$ \tilde{\nabla}^2 \tilde{P}= 0 $$

Which you can descritize into a sparse matrix fairly cleanly or solve by some other methods. All of this is outlined in my thesis, which just happens to be titled "Properties of Stochastic Flow and Permeability of Random Porous Media". I have also implemented a finite-difference solver that solves these equations in 2, 3 and higher dimensions.

The secret horrible hard part is getting all the BC's right :)