begin
    import Pkg as _Pkg
    haskey(ENV, "PLUTO_PROJECT") && _Pkg.activate(ENV["PLUTO_PROJECT"])
    using Revise
    using Test
    using VoronoiFVM
    using GridVisualize
    using CairoMakie
    CairoMakie.activate!(; type = "svg", visible = false)
    GridVisualize.default_plotter!(CairoMakie)
    using SimplexGridFactory, Triangulate
    using ExtendableGrids
    using PlutoUI, HypertextLiteral, UUIDs
end

Outflow boundary conditions

Source

We show how to implment outflow boundary conditions when the velocities at the boundary are calculated by another equation in the system. A typical case is solute transport in porous media where fluid flow is calculated by Darcy's law which defines the convective velocity in the solute transport equation.

Regard the following system of equations in domain $\Omega\subset \mathbb R^d$:

$$\begin{aligned} \nabla \cdot \vec v &=0\\ \vec v&=-k\nabla p\\ \nabla \cdot \vec j &=0\\ \vec j&= - D\nabla c - c\vec v \end{aligned}$$

The variable p can be seen as a the pressure of a fluid in porous medium. c is the concentration of a transported species.

We subdivide the boundary: $\partial\Omega=Γ_{in}\cup Γ_{out}\cup Γ_{noflow}$ abs set

$$\begin{aligned} p=&1 \quad & c&=c_{in} & \text{on}\quad Γ_{in}\\ p=&0 \quad & \vec j\cdot \vec n &= c\vec v \cdot \vec n & \text{on}\quad Γ_{out}\\ \vec v\cdot \vec n &=0 \quad & \vec j\cdot \vec n &= 0 & \text{on}\quad Γ_{noflow}\\ \end{aligned}$$

Discretization data

Base.@kwdef struct FlowTransportData
    k = 1.0
    v_in = 1.0
    c_in = 0.5
    D = 1.0
    Γ_in = 1
    Γ_out = 2
    ip = 1
    ic = 2
end

FlowTransportData

X = 0:0.1:1

0.0:0.1:1.0

darcyvelo(u, data) = data.k * (u[data.ip, 1] - u[data.ip, 2])

darcyvelo (generic function with 1 method)

function flux(y, u, edge, data)
    vh = darcyvelo(u, data)
    y[data.ip] = vh

    bp, bm = fbernoulli_pm(vh / data.D)
    y[data.ic] = data.D * (bm * u[data.ic, 1] - bp * u[data.ic, 2])
end

flux (generic function with 1 method)

function bcondition(y, u, bnode, data)
    boundary_neumann!(y, u, bnode; species = data.ip, region = data.Γ_in, value = data.v_in)
    boundary_dirichlet!(y, u, bnode; species = data.ip, region = data.Γ_out, value = 0)
    boundary_dirichlet!(y, u, bnode; species = data.ic, region = data.Γ_in, value = data.c_in)
end

bcondition (generic function with 1 method)

This function describes the outflow boundary condition. It is called on edges (including interior ones) which have at least one ode on one of the outflow boundaries. Within this function outflownode can be used to identify that node. There is some ambiguity in the case that both nodes are outflow nodes, in that case it is assumed that the contribution is zero. In the present case this is guaranteed by the constant Dirichlet boundary condition for the pressure.

function boutflow(y, u, edge, data)
    y[data.ic] = -darcyvelo(u, data) * u[data.ic, outflownode(edge)]
end

boutflow (generic function with 1 method)

function flowtransportsystem(grid; kwargs...)
    data = FlowTransportData(; kwargs...)
    VoronoiFVM.System(grid;
                      flux,
                      bcondition,
                      boutflow,
                      data,
                      outflowboundaries = [data.Γ_out],
                      species = [1, 2],)
end

flowtransportsystem (generic function with 1 method)

function checkinout(sys, sol)
    data = sys.physics.data
    tfact = TestFunctionFactory(sys)
    tf_in = testfunction(tfact, [data.Γ_out], [data.Γ_in])
    tf_out = testfunction(tfact, [data.Γ_in], [data.Γ_out])
    (; in = integrate(sys, tf_in, sol), out = integrate(sys, tf_out, sol))
end

checkinout (generic function with 1 method)

1D Case

grid = simplexgrid(X)

ExtendableGrids.ExtendableGrid{Float64, Int32};
dim: 1 nodes: 11 cells: 10 bfaces: 2

sys1 = flowtransportsystem(grid);

sol1 = solve(sys1; verbose = "n");

t1 = checkinout(sys1, sol1)

(in = [1.0, 0.5000000000000002], out = [-1.0, -0.5000000000000002])

@test t1.in ≈ -t1.out

Test Passed

@test maximum(sol1[2, :]) ≈ 0.5

Test Passed

@test minimum(sol1[2, :]) ≈ 0.5

Test Passed

2D Case

begin
    g2 = simplexgrid(X, X)
    bfacemask!(g2, [1, 0.3], [1, 0.7], 5)
end

ExtendableGrids.ExtendableGrid{Float64, Int32};
dim: 2 nodes: 121 cells: 200 bfaces: 40

gridplot(g2; size = (300, 300))

sys2 = flowtransportsystem(g2; Γ_in = 4, Γ_out = 5);

sol2 = solve(sys2; verbose = "n")

2×121 Matrix{Float64}:
 1.12521  1.02537  0.925877  0.826948  …  0.453027  0.377299  0.321519  0.298904
 0.5      0.5      0.5       0.5          0.5       0.5       0.5       0.5

let
    vis = GridVisualizer(; size = (700, 300), layout = (1, 2))
    scalarplot!(vis[1, 1], g2, sol2[1, :])
    scalarplot!(vis[1, 2], g2, sol2[2, :]; limits = (0, 1))
    reveal(vis)
end

t2 = checkinout(sys2, sol2)

(in = [1.0000000000000013, 0.5000000000000006], out = [-1.0000000000000007, -0.5000000000000001])

@test t2.in ≈ -t2.out

Test Passed

@test maximum(sol2[2, :]) ≈ 0.5

Test Passed

@test minimum(sol2[2, :]) ≈ 0.5

Test Passed

3D Case

begin
    g3 = simplexgrid(X, X, X)
    bfacemask!(g3, [0.3, 0.3, 0], [0.7, 0.7, 0], 7)
end

ExtendableGrids.ExtendableGrid{Float64, Int32};
dim: 3 nodes: 1331 cells: 6000 bfaces: 1200

gridplot(g3; size = (300, 300))

sys3 = flowtransportsystem(g3; Γ_in = 6, Γ_out = 7);

sol3 = solve(sys3; verbose = "n")

2×1331 Matrix{Float64}:
 0.547438  0.539229  0.516946  0.488884  …  1.32867  1.32914  1.32951  1.32966
 0.5       0.5       0.5       0.5          0.5      0.5      0.5      0.5

let
    vis = GridVisualizer(; size = (700, 300), layout = (1, 2))
    scalarplot!(vis[1, 1], g3, sol3[1, :])
    scalarplot!(vis[1, 2], g3, sol3[2, :]; limits = (0, 1))
    reveal(vis)
end

t3 = checkinout(sys3, sol3)

(in = [1.000000000000037, 0.5000000000000188], out = [-1.0000000000000389, -0.5000000000000194])

@test t3.in ≈ -t3.out

Test Passed

@test maximum(sol3[2, :]) ≈ 0.5

Test Passed

@test minimum(sol3[2, :]) ≈ 0.5

Test Passed

Built with Julia 1.10.2 and

CairoMakie 0.10.12
ExtendableGrids 1.2.3
GridVisualize 1.1.7
HypertextLiteral 0.9.5
Pkg 1.10.0
PlutoUI 0.7.55
Revise 3.5.13
SimplexGridFactory 0.5.20
Triangulate 2.3.2
VoronoiFVM 1.16.0