GitHub

Some problems with Voronoi FVM

Source

Draft. J. Fuhrmann, Oct. 29. 2021. Updated Dec 19, 2021.

We discuss one of the critical cases for application the Voronoi finite volume method. We provide some practical fix and opine that the finite element method probably has the same problems.

# Convert this cell to markdown in order to enable Pluto's inbuilt package managerif isdefined(Main, :PlutoRunner)    using Pkg    docsdir = joinpath(@__DIR__, "..", "docs")    if isdir(docsdir)        Pkg.activate(docsdir)    endend
begin    using Revise    using Test    using VoronoiFVM    using ExtendableGrids    using PlutoUI    using GridVisualize    if isdefined(Main, :PlutoRunner)        using CairoMakie        CairoMakie.activate!(; type = "png")        GridVisualize.default_plotter!(CairoMakie)    endend;

Transient problem

This problem was suggested by R. Eymard.

Regard the following problem coupling Darcy's equation with Fick's law and transport:

$$ \begin{aligned} \vec v &= k \nabla p \\ \nabla \cdot \vec v &= 0\\ \partial_t (\phi c) - \nabla \cdot (D\nabla c + c \vec v) &= 0 \end{aligned}$$

The domain is described by the following discretization grid:

In the center of the domain, we assume a layer with high permeability.

As boundary conditions for the pressure \(p\) we choose fixed pressure values at the left and right boundaries of the domain, triggering a constant pressure gradient throughout the domain.

At the inlet of the high permeability layer, we set \(c=1\), and at the outlet, we set \(c=0\).

The high permeability layer has length L=10 and width W= 0.5.

We solve the time dependent problem on three types of rectangular grids with the same resolution in \(x\) direction and different variants to to handle the high permeability layer.

  • grid_n - a "naive" grid which just resolves the permeability layer and the surrounding material with equally spaced (in y direction) grids

  • grid_1 - a 1D grid of the high permeability layer. With high permeability contrast, the solution of the 2D case at y=0 should coincide with the 1D solution

  • grid_f - a "fixed" 2D grid which resolves the permeability layer and the surrounding material with equally spaced (in y direction) grids and "protection layers" of width ε_fix=0.0001 correcting the size of high permeability control volumes

Results

In the calculations, we ramp up the inlet concentration and measure the amount of solute flowing through the outlet - the breaktrough curve.

nref = 1
1
tend = 100
100
ε_fix = 1.0e-4
0.0001
grid_n, sol_n, bt_n = trsolve(grid_2d(; nref = nref); tend = tend);
sum(bt_n)
18.143158169851787
@test sum(bt_n) ≈ 18.143158169851787
�[32m�[1mTest Passed�[22m�[39m
grid_1, sol_1, bt_1 = trsolve(grid_1d(; nref = nref); tend = tend);
@test sum(bt_1) ≈ 20.66209910195916
�[32m�[1mTest Passed�[22m�[39m
grid_f, sol_f, bt_f = trsolve(grid_2d(; nref = nref, ε_fix = ε_fix); tend = tend);
@test sum(bt_f) ≈ 20.661131375044135
�[32m�[1mTest Passed�[22m�[39m
grid_ϕ, sol_ϕ, bt_ϕ = trsolve(grid_2d(; nref = nref); ϕ = [1.0e-3, 1], tend = tend);
@test sum(bt_ϕ) ≈ 20.412256299447236
�[32m�[1mTest Passed�[22m�[39m
begin    p1 = GridVisualizer(;        resolution = (500, 200),        xlabel = "t",        ylabel = "outflow",        legend = :rb,        title = "Breakthrough Curves"    )    scalarplot!(p1, sol_n.t, bt_n; label = "naive grid", color = :red)    scalarplot!(        p1,        sol_1.t,        bt_1;        label = "1D grid",        markershape = :x,        markersize = 10,        clear = false,        color = :green    )    scalarplot!(        p1,        sol_f.t,        bt_f;        label = "grid with fix",        markershape = :circle,        color = :green,        clear = false    )    scalarplot!(        p1,        sol_ϕ.t,        bt_ϕ;        label = "modified ϕ",        markershape = :cross,        color = :blue,        clear = false    )    reveal(p1)end

Here, we plot the solutions for the grid_n case and the grid_f case.

Time: 10.0

scalarplot(grid_n, sol_n(t)[ic, :]; resolution = (500, 200), show = true)
scalarplot(grid_f, sol_f(t)[ic, :]; resolution = (500, 200), show = true)

Reaction-Diffusion problem

Here we solve the following problem:

$$ -\nabla D \nabla u + R u = 0$$

where D is large in the high permeability region and small otherwise. R is a constant.

Results

rdgrid_1, rdsol_1, of_1 = rdsolve(grid_1d(; nref = nref));
@test of_1 ≈ -0.013495959676585267
�[32m�[1mTest Passed�[22m�[39m
rdgrid_n, rdsol_n, of_n = rdsolve(grid_2d(; nref = nref));
@test of_n ≈ -0.00023622450350365264
�[32m�[1mTest Passed�[22m�[39m
rdgrid_f, rdsol_f, of_f = rdsolve(grid_2d(; nref = nref, ε_fix = ε_fix));
@test of_f ≈ -0.013466874615165499
�[32m�[1mTest Passed�[22m�[39m
rdgrid_r, rdsol_r, of_r = rdsolve(grid_2d(; nref = nref); R = [0, 0.1]);
@test of_r ≈ -0.013495959676764535
�[32m�[1mTest Passed�[22m�[39m

We measure the outflow at the outlet. As a result, we obtain:

  • 1D case: -0.013495959676585248

  • 2D case, naive grid: -0.000236224503503653

  • 2D case, grid with "protective layer": -0.013466874615165499

  • 2D case, naive grid, "modified" R: -0.013495959676764509

scalarplot(rdgrid_1, rdsol_1; resolution = (300, 200))
scalarplot(rdgrid_n, rdsol_n; resolution = (500, 200))
scalarplot(rdgrid_f, rdsol_f; resolution = (500, 200))

Discussion

Transient case

As there will be nearly no flow in y-direction, we should get the very same results in all four cases for small permeability values in the low permeability region.

In the grid_n case, the heterogeneous control volumina ovrestimate the storage capacity which shows itself in a underestimation of the transferred solute.

With the high permeability contrast, the results for heterogeneous domain should be essentially equal to those for 1D domain. However, with a coarse resolution in y-direction, we see large differences in the transient behaviour of the breaktrough curve compared to the 1D case. The introduction of a thin protection layer leads to reasonable results.

Alternatively, the porosity of the low permeability region can be modified. Arguably, this is the case in practice, see e.g. Ackerer et al, Transport in Porous Media35:345–373, 1999 (eq. 7).

Reaction diffusion case

In this case, we look at a homogeneous reaction in a domain divided into a high and low diffusion region. With high contrast in the diffusion coefficients, the reasonable assumption is that the reaction takes place only in the high diffusion region, and the un-consumed share of species leaves at the outlet.

In this case we observe a similar related problem which can be fixed by adding a thin layer of control volumes at the boundary. No problem occurs if the reaction rate at in the low diffusion region is zero.

Conclusion

Here, we indeed observe problem with the Voronoi approach: care must be taken to handle the case of hetero interfaces in connection with transient processes and/or homogeneous reactions. In these cases it should be analyzed if the problem occurs, and why, and it appears, that the discussion should not be had without reference to the correct physical models. A remedy based on meshing is available at least for straight interfaces.

Opinion

With standard ways of using finite elements, the issue described here will occur in a similar way, so the discussion is indeed with the alternative "cell centered" finite volume approach which places interfaces at the boundaries of the control volumes rather than along the edges of a underlying triangulation.

Drawbacks of two point flux Voronoi methods based on simplicial meshes (as tested here):

  • Anisotropic diffusion is only correct with aligned meshes

  • Reliance on boundary conforming Delaunay property of the underlying mesh, thus narrowing the available meshing strategies

  • The issue described in the present notebook. However, in both cases discussed here, IMHO it might "go away" depending on the correct physics. There should be more discussions with relevant application problems at hand.

Advantages (compared to the cell centered approach placing collocation points away from interfaces)

  • Availability of P1 interpolant on simplices for visualization, interpolation, coupling etc.

  • Mesh generators tend to place interfaces at triangle edges.

  • Dirichlet BC can be applied exactly

  • There is a straightforward way to link interface processes with bulk processes, e.g. an adsorption reaction is easily described by a reaction term at the boundary which involves interface and bulk value available at the same mesh node.

Appendix

Domain data

Sizes:

begin    L = 10   # length of the high perm layer    W = 0.5  # width of high perm layer    Wlow = 2 # width of adjacent low perm layersend;

Boundary conditions:

begin    const Γ_top = 3    const Γ_bot = 1    const Γ_left = 4    const Γ_right = 2    const Γ_in = 5    const Γ_out = 2end;
begin    Ω_low = 1    Ω_high = 2end;
function grid_2d(; nref = 0, ε_fix = 0.0)    nx = 10 * 2^nref    ny = 1 * 2^nref    nylow = 3 * 2^nref    xc = linspace(0, L, nx + 1)    y0 = linspace(-W / 2, W / 2, ny + 1)    if ε_fix > 0.0        yfix = [W / 2, W / 2 + ε_fix]        ytop = glue(yfix, linspace(yfix[end], Wlow, nylow + 1))    else        ytop = linspace(W / 2, Wlow, nylow + 1)    end    yc = glue(-reverse(ytop), glue(y0, ytop))    grid = simplexgrid(xc, yc)    cellmask!(grid, [0, -W / 2], [L, W / 2], Ω_high)    bfacemask!(grid, [0, -W / 2], [0, W / 2], Γ_in)    return bfacemask!(grid, [L, -W / 2], [L, W / 2], Γ_out)end
grid_2d (generic function with 1 method)
function grid_1d(; nref = 0)    nx = 10 * 2^nref    xc = linspace(0, L, nx + 1)    grid = simplexgrid(xc)    cellmask!(grid, [0], [L], Ω_high)    bfacemask!(grid, [0], [0], Γ_in)    bfacemask!(grid, [L], [L], Γ_out)    return gridend
grid_1d (generic function with 1 method)

Transient solver

Pressure index in solution

const ip = 1;

Concentration index in solution

const ic = 2;

Generate breaktrough courve from transient solution

function breakthrough(sys, tf, sol)    of = similar(sol.t)    t = sol.t    of[1] = 0    for i in 2:length(sol.t)        of[i] = -integrate(sys, tf, sol.u[i], sol.u[i - 1], t[i] - t[i - 1])[ic]    end    return ofend
breakthrough (generic function with 1 method)

Transient solver:

function trsolve(        grid;        κ = [1.0e-3, 5],        D = [1.0e-12, 1.0e-12],        Δp = 1.0,        ϕ = [1, 1],        tend = 100,    )    function flux(y, u, edge, data)        y[ip] = κ[edge.region] * (u[ip, 1] - u[ip, 2])        bp, bm = fbernoulli_pm(y[ip] / D[edge.region])        y[ic] = D[edge.region] * (bm * u[ic, 1] - bp * u[ic, 2])        return nothing    end    function stor(y, u, node, data)        y[ip] = 0        y[ic] = ϕ[node.region] * u[ic]        return nothing    end    dim = dim_space(grid)    function bc(y, u, bnode, data)        c0 = ramp(bnode.time; dt = (0, 0.001), du = (0, 1))        boundary_dirichlet!(y, u, bnode, ic, Γ_in, c0)        boundary_dirichlet!(y, u, bnode, ic, Γ_out, 0)        boundary_dirichlet!(y, u, bnode, ip, Γ_in, Δp)        boundary_dirichlet!(y, u, bnode, ip, Γ_out, 0)        if dim > 1            boundary_dirichlet!(y, u, bnode, ip, Γ_left, Δp)            boundary_dirichlet!(y, u, bnode, ip, Γ_right, 0)        end        return nothing    end    sys = VoronoiFVM.System(        grid;        check_allocs = true,        flux = flux,        storage = stor,        bcondition = bc,        species = [ip, ic],    )    inival = solve(sys; inival = 0, time = 0.0)    factory = TestFunctionFactory(sys)    tfc = testfunction(factory, [Γ_in, Γ_left, Γ_top, Γ_bot], [Γ_out])    sol = VoronoiFVM.solve(        sys;        inival = inival,        times = [0, tend],        Δt = 1.0e-4,        Δt_min = 1.0e-6,    )    bt = breakthrough(sys, tfc, sol)    if dim == 1        bt = bt * W    end    return grid, sol, btend
trsolve (generic function with 1 method)

Reaction-Diffusion solver

function rdsolve(grid; D = [1.0e-12, 1.0], R = [1, 0.1])    function flux(y, u, edge, data)        y[1] = D[edge.region] * (u[1, 1] - u[1, 2])        return nothing    end    function rea(y, u, node, data)        y[1] = R[node.region] * u[1]        return nothing    end    function bc(y, u, bnode, data)        boundary_dirichlet!(y, u, bnode, 1, Γ_in, 1)        boundary_dirichlet!(y, u, bnode, 1, Γ_out, 0)        return nothing    end    sys = VoronoiFVM.System(        grid;        flux = flux,        reaction = rea,        species = 1,        bcondition = bc,        check_allocs = true    )    dim = dim_space(grid)    sol = VoronoiFVM.solve(sys)    factory = TestFunctionFactory(sys)    tf = testfunction(factory, [Γ_in, Γ_left, Γ_top, Γ_bot], [Γ_out])    of = integrate(sys, tf, sol)    fac = 1.0    if dim == 1        fac = W    end    return grid, sol[1, :], of[1] * facend
rdsolve (generic function with 1 method)