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 articifial diffusion.

The computed solution for the default parameters looks like this:

module Example265_FlowTransport

using ExtendableFEM
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
end
function c_inlet!(result, qpinfo)
    result[1] = (1-qpinfo.x[2])*qpinfo.x[2]
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])
end

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

function kernel_inlet!(result, input, qpinfo)
    c_inlet!(result, qpinfo)
    result[1] *= -input[1]
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)
        τ = 1e3
        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 = 1e-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 = 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 = 1e-12, constant_matrix = true, kwargs...)

    # print minimal and maximal concentration to check max principle (shoule 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)

    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 = 1 : num_cells(xgrid)
            nfaces4cell = num_targets(cellfaces, cell)
            for cf = 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,1e-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)
    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.