I have a small 2D system I'm trying to model using a non-linear extension of Darcy's law for fluid flow in porous media. I'm primarily interested in the local flow velocity, not necessarily the pressure, that's later used in the real model I'm working on:
The original equation known as Darcy's Law used in ground water modeling is essentially Poisson's equation, where $\vec{q}$ is the fluid velocity, $\mu$ and $\kappa$ are constants that together describe the hydraulic conductivity of the porous medium, $p$ is local hydrostatic head and f is a source:
$$ -\frac{\mu}{\kappa} \nabla p = \vec{q} $$
$$ \nabla \cdot ( \vec{q} ) = f = -\frac{\mu}{\kappa} \nabla^2 p $$
The non-linear equation I want to solve instead is the following, known as the Darcy–Forchheimer law. As you can see it's essentially a polynomial instead of a linear model like Darcy's law:
$$ \frac {\partial p}{\partial x} = \frac {\mu }{\kappa} \vec{q} - \frac {\rho }{\kappa _{1}}\vec{q}^{2} $$
I've picked up on some authors [1] that go through the calculus and have figured out a way to represent $\vec{q}$ as a function of $p$ by calculating the inverse of the flow/pressure equation, and I see somewhat how they arrive at their methodology, but what I'm interested in is actually calculating the flux ($\vec{q}$) itself, which I use elsewhere in the model I'm developing - I don't have much of a use for $p$ besides in the initial problem setup with some simple Dirichlet boundary conditions.
Somewhat related, I've also read that for some numerical methods such as DG FEM, diffusion dominated problems can very unstable due to the non-directional flux. In this paper [2], the authors investigate some 'equation splitting' methods where they solve one equation for the state variable, and another for the flux, essentially the first two equations I showed above. This apparently relieves some of the instability issues caused by diffusion
Looking at these two, it looks like I have two options:
- Use the methodology of [1] to calculate the pressure throughout the system, and then calculate the gradient of the hydrostatic head to generate the local velocity, Or...
- Split my equation into two equations and solve them together. I end up having the solve two equations at once, but I need to have the velocity regardless so that doesn't seem so bad. Other than that I know nothing about this method.
I don't know how to implement the equation splitting method as I've never done it before, so I have two questions:
- Is this even a good idea? Is splitting the equations up into double the number of equations disadvantageous for any other reason than doubling the number of state equations to solve?
- Is there a source that goes into the actual implementation of equation splitting for (preferably) the finite element or finite difference methods. I am thoroughly intimidated by any DG theory I come across, and I'm hoping there is an easier way to solve this problem
[1] https://arxiv.org/pdf/1508.00294.pdf
[2] https://www.cs.utah.edu/~kirby/Publications/Kirby-11.pdf
