Introduction to FEniCS - Darcy

Introduction to FEniCS - Darcy#

../_images/fenics_logo.png

In this notebook we will introduce the concept of pressure driven flow in porous media, as described by Darcy’s equation. You will learn how to:

  • Solve for a Darcian pressure field

  • Post-process a velocity field for application to fluid flow in porous media

Porous flow and Darcy’s equation in steady state#

../_images/Darcy.png

Mathematical formulation:#

Darcy’s law is an equation that describes the flow of a fluid through a porous medium. In the absence of gravitational forces, it provides a simple relationship between the flow rate through a porous medium and the pressure drop over a given distance. For homogeneous permeable medium with permeability \({\kappa}\) and fluid dynamic viscosity \({\mu}\), one we can derive:

\[ q = - \: \frac {\kappa} {\mu} \: \nabla p \:\:\:\:\:\:\:\:\:\:\:\: (1) \]

where \(q\) is the fluid discharge per unit area (m/s) and \(\nabla p\) is the fluid pressure gradient vector (Pa/m). The discharge flux is often referred to as the Darcy flux or Darcy velocity, and is an average quantity of the fluid front through the porous media, it is different from the local fluid velocity traveling through the pore network. The fluid velocity \(v_f\) can be recovered by rescaling the Darcy flux \(q\) by the porosity of the media \(\phi\).

\[ v_f = \frac {q} {\phi} \]

As now, we have one equation for two unknowns (the Darcy flux and the pressure gradient). By using the mass balance equation of the fluid, one can derive an expression for the fluid pressure:

\[ \frac {\partial \: \phi\rho_f} {\partial t} = \nabla \cdot \rho_f v_f \]

Which simplifies for constant fluid density \(\rho_f\) and constant porosity field \(\phi\) to:

\[ \nabla \cdot v_f = 0 \:\:\:\:\:\: or \:\:\:\:\:\: \nabla \cdot q = 0 \:\:\:\:\:\:\:\:\:\:\:\: (2) \]

Using equation (1) in (2) gives:

\[ \nabla \cdot ( - \: \frac {\kappa} {\mu} \: \nabla p ) = 0 \]

We have derived a Poisson equation as an expression for the fluid pressure field. One can solve this equation given some boundary conditions.

Implementation#

We want to implement this variational problem

%matplotlib inline

# Import all the functions from DOLFIN (computational backend of the FEniCS project)
from dolfin import *

Discretizing the domain and space of functions#

We can now use FEniCS to solve our problem. First, we discretize our domain by creating a mesh of \(\Omega\).

# 20 is number of intervals Omega is divided intoc
mesh = UnitSquareMesh(20, 20)
plot(mesh)

elem = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, elem)
../_images/b7b49a0f1f4ed483467151814d7da7b112dfbd8f8cb4601e6898ce0100f0246a.png

Specifying the weak formulation#

Next, all we need to do is tell FEniCS the weak formulation of our problem.

kappa_over_mu = Constant(1.0)  # physical material property
phi = Constant(0.1)  # porosity, ranging from 0 to 1
S = Constant(0.0)  # source term
v = TestFunction(W)
p = TrialFunction(W)

Set the boundary conditions#

Next, all we need to do is tell FEniCS the weak formulation of our problem.

value_left = Constant(1.0)
value_right = Constant(0.0)

# test with differente boundary conditions...
# Imposing Dirichlet BC to the left boundary node
bc_l = DirichletBC(W, value_left, "on_boundary && near(x[0], 0)")
# Imposing Dirichlet BC to the right boundary node
bc_r = DirichletBC(W, value_right, "on_boundary && near(x[0], 1)")
bcs = [bc_l, bc_r]   # list of boundary conditions to apply to the problem

Exercise 1: Create and solve the Darcy equation#

# Add solution here

Solution:#

Hide code cell source 
F = kappa_over_mu * dot(grad(p), grad(v)) * dx - S * \
    v * dx  # residual form of our equation
a, L = system(F)
ph = Function(W)  # place to store the solution
solve(a == L, ph, bcs)
plot(ph)
Hide code cell output 
Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Solving linear variational problem.

<matplotlib.tri._tricontour.TriContourSet at 0x7f828b697f10>
../_images/67366981f76fd850657318faf3aa3dfb75dd95a5db7aab504740f6429f796d75.png

Estimate the velocity field from the pressure field (postprocessing)#

In the previous tutorials, we have only worked with scalar quantities (such as concentration, pressure, temperature). However, as we would like to get the velocity field $\( v_f = \frac {q} {\phi} \)$ we need to define a suitable finite element space for this quantity.

As \(q\) is a first order Lagrange function, it’s derivative is a cell-wise continuous function (as explained in Lecture 1). Thus, we create a vector-element, which will have the same number of dimensions as the mesh (1 in 1D, 2 in 2D, 3 in 3D).

# This is new type of element (for vector field discretization)
elem_v = VectorElement("DG", mesh.ufl_cell(), 0)
W_v = FunctionSpace(mesh, elem_v)

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Next, we need to compute \(v_f\). We do this by a projection, which means solving the linear problem: Find \(v_f\in W_v\) such that $\( \int_\Omega (v_f - \frac{q}{\phi})v~\mathrm{d}x = 0\qquad\forall v\in W_v \)$ There are several ways of computing this projection. The simplest (but most expensive if we need to do it multiple times) is to use the project function in FEniCS.

import matplotlib.pyplot as plt
vf = project(-kappa_over_mu * grad(ph) / phi, W_v)
fig = plt.figure()
im = plot(vf)
plt.colorbar(im, format="%.2e")

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

<matplotlib.colorbar.Colorbar at 0x7f8282086ad0>
../_images/dff7ebffba13593e9fddf2caf2c59d7429622845de625875ff23e34c38e0a0e9.png

Paraview#

We can output the solutions to Paraview for post-processing. We do this by using the xdmf-format, as shown below

with XDMFFile("Darcy_solutions.xdmf") as xdmf:
    xdmf.write(vf, 0.0)
    xdmf.write(ph, 0.0)

Advanced: Reusable projection#

As mentioned above, a projection is solve the problem

\[ a(u, v) = l(v) \]

where

\[ a(u,v)=\int_\Omega u\cdot v~\mathrm{d}x \]
\[ l(v)=\int_\Omega f\cdot v~\mathrm{d}x \]
import matplotlib.pyplot as plt


class Projector():
    def __init__(self, V):
        u = TrialFunction(V)
        v = TestFunction(V)
        a = inner(u, v) * dx
        self.A = assemble(a)
        self.u = Function(V)
        self.solver = LUSolver(self.A)

    def __call__(self, f):
        v = TestFunction(self.u.function_space())
        b = assemble(inner(f, v)*dx)
        self.solver.solve(self.u.vector(), b)
        return self.u


projector = Projector(W_v)
vel_f = projector(-kappa_over_mu * grad(ph) / phi)
fig = plt.figure()
im = plot(vel_f)
plt.colorbar(im, format="%.2e")

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

<matplotlib.colorbar.Colorbar at 0x7f8281f5f430>
../_images/dff7ebffba13593e9fddf2caf2c59d7429622845de625875ff23e34c38e0a0e9.png

Advanced: Local projection#

As we are using a discontinuous function space to represent the velocity the degrees of freedom in each cell is independent of the others. Therefore, we can create a projection function using a LocalSolver.

def local_project(f, V):
    u = TrialFunction(V)
    v = TestFunction(V)
    a_proj = inner(u, v)*dx
    b_proj = inner(f, v)*dx
    solver = LocalSolver(a_proj, b_proj)
    solver.factorize()
    u = Function(V)
    solver.solve_local_rhs(u)
    return u


vel_f = local_project(-kappa_over_mu * grad(ph) / phi, W_v)
fig = plt.figure()
plot(mesh, linewidth=0.1)
im = plot(vel_f)
plt.colorbar(im, format="%.2e")

<matplotlib.colorbar.Colorbar at 0x7f8280e62f80>
../_images/b0e50786b55c0c42bd057c091ac5f785f4bee6286e1c7c1993638108faea8e8c.png
Contents