265 : Flow + Transport

(source code)

This example solve the Stokes problem in an Omega-shaped pipe and then uses the velocity in a transport equation for a species with a certain inlet concentration. Altogether, we are looking for a velocity $\mathbf{u}$, a pressure $\mathbf{p}$ and a stationary species concentration $\mathbf{c}$ such that

\[\begin{aligned} - \mu \Delta \mathbf{u} + \nabla p & = 0\\ \mathrm{div}(\mathbf{u}) & = 0\\ \mathbf{c}_t - \kappa \Delta \mathbf{c} + \mathbf{u} \cdot \nabla \mathbf{c} & = 0 \end{aligned}\]

with some viscosity parameter and diffusion parameter $\kappa$.

The diffusion coefficient for the species is chosen (almost) zero such that the isolines of the concentration should stay parallel from inlet to outlet. For the discretisation of the convection term in the transport equation three possibilities can be chosen:

  1. Classical Bernardi–Raugel stationary finite element discretisations $\mathbf{u}_h \cdot \nabla \mathbf{c}_h$ [set FVtransport = false, reconstruct = false]
  2. As in 1. but with divergence-free reconstruction operator in convection term $\Pi_\text{reconst} \mathbf{u}_h \cdot \nabla \mathbf{c}_h$ [set FVtransport = false, reconstruct = true]
  3. Time-dependent upwind finite volume discretisation for $\kappa = 0$ based on normal fluxes along the faces [set FVtransport = true]

Observe that the divergence-free postprocessing helps a lot for mass conservation, but is still not perfect. The finite volume upwind discretisation ensures mass conservation.

Note, that the transport equation is very convection-dominated and no stabilisation in the finite element discretisations was used here (but instead a nonzero $\kappa$). Also note, that only the finite volume discretisation perfectly obeys the maximum principle for the concentration but the isolines do no stay parallel until the outlet is reached, possibly due to artificial diffusion.

The computed solution for the default parameters looks like this:

module Example265_FlowTransport

using ExtendableFEM
using ExtendableFEMBase
using ExtendableGrids
using SimplexGridFactory
using Triangulate

# boundary data
function u_inlet!(result, qpinfo)
    x = qpinfo.x
    result[1] = 4 * x[2] * (1 - x[2])
    result[2] = 0
    return nothing
end
function c_inlet!(result, qpinfo)
    result[1] = (1 - qpinfo.x[2]) * qpinfo.x[2]
    return nothing
end

function kernel_stokes_standard!(result, u_ops, qpinfo)
    ∇u, p = view(u_ops, 1:4), view(u_ops, 5)
    μ = qpinfo.params[1]
    result[1] = μ * ∇u[1] - p[1]
    result[2] = μ * ∇u[2]
    result[3] = μ * ∇u[3]
    result[4] = μ * ∇u[4] - p[1]
    result[5] = -(∇u[1] + ∇u[4])
    return nothing
end

function kernel_convection!(result, ∇T, u, qpinfo)
    result[1] = ∇T[1] * u[1] + ∇T[2] * u[2]
    return nothing
end

function kernel_inlet!(result, input, qpinfo)
    c_inlet!(result, qpinfo)
    result[1] *= -input[1]
    return nothing
end


# everything is wrapped in a main function
function main(; nrefs = 4, Plotter = nothing, reconstruct = true, FVtransport = true, parallel = false, npart = 8, μ = 1, kwargs...)

    # load mesh and refine
    xgrid = uniform_refine(
        simplexgrid(
            Triangulate;
            points = [0 0; 3 0; 3 -3; 7 -3; 7 0; 10 0; 10 1; 6 1; 6 -2; 4 -2; 4 1; 0 1]',
            bfaces = [1 2; 2 3; 3 4; 4 5; 5 6; 6 7; 7 8; 8 9; 9 10; 10 11; 11 12; 12 1]',
            bfaceregions = [1; 1; 1; 1; 1; 2; 3; 3; 3; 3; 3; 4],
            regionpoints = [0.5 0.5;]',
            regionnumbers = [1],
            regionvolumes = [1.0]
        ), nrefs
    )

    if parallel
        xgrid = partition(xgrid, RecursiveMetisPartitioning(npart = npart))
    end

    # define unknowns
    u = Unknown("u"; name = "velocity", dim = 2)
    p = Unknown("p"; name = "pressure", dim = 1)
    T = Unknown("T"; name = "temperature", dim = 1)

    id_u = reconstruct ? apply(u, Reconstruct{HDIVBDM1{2}, Identity}) : id(u)

    # define first sub-problem: Stokes equations to solve for velocity u
    PD = ProblemDescription("Stokes problem")
    assign_unknown!(PD, u)
    assign_unknown!(PD, p)
    assign_operator!(PD, BilinearOperator(kernel_stokes_standard!, [grad(u), id(p)]; params = [μ], parallel = parallel, kwargs...))
    assign_operator!(PD, InterpolateBoundaryData(u, u_inlet!; regions = 4, parallel = parallel, kwargs...))
    assign_operator!(PD, HomogeneousBoundaryData(u; regions = [1, 3], kwargs...))

    # add transport equation of species
    PDT = ProblemDescription("transport problem")
    assign_unknown!(PDT, T)
    if FVtransport ## FVM discretisation of transport equation (pure upwind convection)
        τ = 1.0e3
        assign_operator!(PDT, CallbackOperator(assemble_fv_operator!(), [u]; kwargs...))
        assign_operator!(PDT, BilinearOperator([id(T)]; store = true, factor = 1 / τ, parallel = parallel, kwargs...))
        assign_operator!(PDT, LinearOperator([id(T)], [id(T)]; factor = 1 / τ, parallel = parallel, kwargs...))
    else ## FEM discretisation of transport equation (with small diffusion term)
        assign_operator!(PDT, BilinearOperator([grad(T)]; factor = 1.0e-6, parallel = parallel, kwargs...))
        assign_operator!(PDT, BilinearOperator(kernel_convection!, [id(T)], [grad(T)], [id_u]; parallel = parallel, kwargs...))
        assign_operator!(PDT, InterpolateBoundaryData(T, c_inlet!; regions = [4], kwargs...))
    end

    # generate FESpaces and a solution vector for all 3 unknowns
    FETypes = [H1BR{2}, L2P0{1}, FVtransport ? L2P0{1} : H1P1{1}]
    FES = [FESpace{FETypes[j]}(xgrid) for j in 1:3]
    sol = FEVector(FES; tags = [u, p, T])

    # solve the two problems separately
    sol = solve(PD; init = sol, kwargs...)
    sol = solve(PDT; init = sol, maxiterations = 20, target_residual = 1.0e-12, constant_matrix = true, kwargs...)

    # print minimal and maximal concentration to check max principle (should be in [0,1])
    println("\n[min(c),max(c)] = [$(minimum(view(sol[T]))),$(maximum(view(sol[T])))]")

    # plot
    plt = plot([id(u), id(T)], sol; Plotter = Plotter, ncols = 1, rasterpoints = 40, width = 800, height = 800)

    return sol, plt
end

# pure convection finite volume operator for transport
function assemble_fv_operator!()

    BndFluxIntegrator = ItemIntegrator(kernel_inflow!, [normalflux(1)]; entities = ON_BFACES)
    FluxIntegrator = ItemIntegrator([normalflux(1)]; entities = ON_FACES)
    fluxes::Matrix{Float64} = zeros(Float64, 1, 0)

    return function closure(A, b, args; assemble_matrix = true, assemble_rhs = true, kwargs...)

        # prepare grid and stash
        xgrid = args[1].FES.xgrid
        nfaces = size(xgrid[FaceCells], 2)
        if size(fluxes, 2) < nfaces
            fluxes = zeros(Float64, 1, nfaces)
        end

        # right-hand side = boundary inflow fluxes if velocity points inward
        if assemble_rhs
            fill!(fluxes, 0)
            evaluate!(fluxes, BndFluxIntegrator, [args[1]])
            facecells = xgrid[FaceCells]
            bface2face = xgrid[BFaceFaces]
            for bface in 1:lastindex(bface2face)
                b[facecells[1, bface2face[bface]]] -= fluxes[bface]
            end
        end

        # assemble upwind finite volume fluxes over cell faces into matrix
        if assemble_matrix
            # integrate normalfux of velocity
            fill!(fluxes, 0)
            evaluate!(fluxes, FluxIntegrator, [args[1]])

            cellfaces = xgrid[CellFaces]
            cellfacesigns = xgrid[CellFaceSigns]
            for cell in 1:num_cells(xgrid)
                nfaces4cell = num_targets(cellfaces, cell)
                for cf in 1:nfaces4cell
                    face = cellfaces[cf, cell]
                    other_cell = facecells[1, face]
                    if other_cell == cell
                        other_cell = facecells[2, face]
                    end
                    flux = fluxes[face] * cellfacesigns[cf, cell]
                    if (other_cell > 0)
                        flux *= 1 // 2 # because it will be accumulated on two cells
                    end
                    if flux > 0 # flow from cell to other_cell or out of domain
                        _addnz(A, cell, cell, flux, 1)
                        if other_cell > 0
                            _addnz(A, other_cell, cell, -flux, 1)
                            # otherwise flow goes out of domain
                        end
                    else # flow from other_cell into cell or into domain
                        _addnz(A, cell, cell, 1.0e-16, 1) # add zero to keep pattern for LU
                        if other_cell > 0 # flow comes from neighbour cell
                            _addnz(A, other_cell, other_cell, -flux, 1)
                            _addnz(A, cell, other_cell, flux, 1)
                        end
                        # otherwise flow comes from outside into domain, handled in rhs side loop above
                    end
                end
            end
        end
        return nothing
    end
end


function kernel_inflow!(result, input, qpinfo)
    return if input[1] < 0 # if velocity points into domain
        c_inlet!(result, qpinfo)
        result[1] *= input[1]
    else
        result[1] = 0
    end
end

end # module

This page was generated using Literate.jl.