# Finite difference scheme for compressible nonisothermal flow in porous media

My challenge is to solve the following system of equations, which describe gas combustion in porous media:

1) Continuity

$\varepsilon \frac{\partial \rho_g}{\partial t} +\frac{\partial}{\partial x} \left(\rho_g u_x\right)=0$

2) Darcy law (momentum)

$u_x=-\frac{k}{\mu} \frac{\partial p}{\partial x}$

3) Equation of state, note the variable temperature

$\rho_g=\frac{M_Rp}{RT_g(x)}$

4) Energy equation for the gas.

5) Energy equation for the solid phase

I have succesfully descritised and solved the case where the velocity, pressure and density are assumed constant, i.e the first three equations drop out. But solving the gasodynamical part proved to be a problem.

Applying an upwind scheme to 1) (as was suggested here: A good finite difference for the continuity equation) yeilds a really harsh stability criteria on the timestep, I am forced to have it as low as 1e-6 with a 1e-2 spacial timestep, even when I take the isothermal case, disregarding the combustion for the time being. And I need at least 1e-3 to resolve the energy equations.

The first three equations can be also coupled together to form

6) $\frac{\partial p}{\partial t} +C\frac{\partial^2}{\partial x^2} \left(p^2 \right)=0$

but only in the isothermal case, so that is of little help.

I know that people have solved 1)-5) and 6) before, but I couldn't find a description of the schemes they used. I tryed searching articles on compressible flow in porous media specifically, but those all deal with way more complex models (multiphase, deformable solids etc) and use very complicated solving methods.

Could someone suggest a good FD scheme for (1)-(3) or say how the stability criteria is formed if one just uses upwinding like I did?

• Is there a reason why you absolutely need to use finite difference method? Fluid dynamics simulations are more naturally resolved using finite volume method because it is naturally conservative.
– Paul
May 18, 2012 at 4:55
• Well, yes, the reason is that changing the method would mean rewriting all the code I have now, and I have a pretty strict deadline on this project, less then a week. I'll be able to report things anyway, but I want to get to the bottom of things :). Feel free to post a finite volume solution though!
– tiam
May 18, 2012 at 10:55
• @Paul True, but only if OP is working on a non-uniform grid. In case of a uniform rectangular grid, finite volume discretization degenerates to finite differences. My opinion is that, if application permitting, FD are great for learning basics, and then FV are the next step. May 20, 2012 at 1:56
• @Paul/@IRO-bot It's more subtle than that. High order conservative finite difference methods exist. There are so many equivalences, especially for simple methods, that to some extent, the choice is only meaningful once we ask what components of the method you would like to stay fixed as you extend the method in a given direction. May 20, 2012 at 5:05
• Usually the equation of state relates pressure, temperature, and density. But I don't see pressure anywhere. And what is $M_{RP}$? May 20, 2012 at 18:43

For hyperbolic problems solved explicitly, you must satisfy a CFL (Courant, Friedrichs, Lewy) condition. This ensures that the scheme will use data only from the domain of dependence for the differential equation. For more information on the CFL condition and why it is necessary, you can read e.g. pp.215-218 in Finite Difference Methods for Ordinary and Partial Differential Equations by LeVeque. Also in that chapter, there are some alternative methods presented but the stability criteria is still there. The CFL condition for the upwind method in 1D is $\frac{u\Delta t}{\Delta x} < 1$.

In Discontinuous Galerkin (DG) methods the use of strong-stability preserving (SSP) Runge-Kutta methods is becoming popular. Though these are explicit methods, they often allow some multiple of the usual Courant number to be used unlike a simple forward Euler method. That means you can take longer time steps but at a greater cost per time step. It might be possible to adapt SSPRK methods to your problem, but I have only seen it done for DG methods and my understanding of their applicability is limited.

It might be possible to use an implicit in time method since they are unconditionally stable. In order to keep the accuracy at an acceptable level, you may end up back with the original restriction on the time step. LeVeque's book seems to suggest that using backward Euler or an Adams method for the time discretization and central differencing for the spatial derivative might work.

I second the vote for Finite Volume methods or, if you want a challenge, DG methods.