260 : Axisymmetric Stokes

(source code)

This example solves the 3D stagnation point flow via the 2.5D axisymmetric formulation of the Navier–Stokes problem that seeks a velocity $\mathbf{u} = (u_z, u_r)$ and pressure $p$ such that

\[\begin{aligned} - \mu\left(\partial^2_r + r^{-1} \partial_r + \partial^2_z - r^{-2} \right) u_r + (u_r \partial_r + u_z \partial_z) u_r + \partial_r p & = \mathbf{f}_r\\ - \mu\left(\partial^2_r + r^{-1} \partial_r + \partial^2_z \right) u_z + (u_r \partial_r + u_z \partial_z) u_z + \partial_z p & = \mathbf{f}_z\\ (\partial_r + r^{-1})u_r + \partial_z u_z & = 0 \end{aligned}\]

with exterior force $\mathbf{f}$ and some viscosity parameter $\mu$.

The axisymmetric formulation assumes that the velocity in some 3D-domain, that is obtained by rotation of a 2D domain $\Omega$, only depends on the distance $r$ to the rotation axis and the $z$-coordinate tangential to the x-axis, but not on the angular coordinate of the cylindric coordinates. The implementation employs $r$-dependent bilinear forms and a Cartesian grid for the 2D $(z,r)$ domain that is assumed to be rotated around the $r=0$-axis.

This leads to the weak formulation

\[\begin{aligned} a(u,v) + b(p,v) & = (f,v) \\ b(q,u) & = 0 \end{aligned}\]

with the bilinear forms

\[\begin{aligned} a(u,v) := \int_{\Omega} \left( \nabla u : \nabla v + r^{-2} u_r v_r \right) r dr dz\\ b(q,v) := \int_{\Omega} q \left( \mathrm{div}(v) + r^{-1} u_r \right) r dr dz \end{aligned}\]

where the usual Cartesian differential operators can be used. The factor $2\pi$ from the integral over the rotation angle drops out on both sides.

The computed solution for the default parameters looks like this:

module Example260_AxisymmetricNavierStokesProblem

using ExtendableFEM
using ExtendableGrids
using SimplexGridFactory
using Triangulate


function kernel_convection!(result, input, qpinfo)
    u, ∇u = view(input, 1:2), view(input, 3:6)
    r = qpinfo.x[1]
    result[1] = r * (∇u[1] * u[1] + ∇u[2] * u[2])
    result[2] = r * (∇u[3] * u[1] + ∇u[4] * u[2])
    return nothing
end

function kernel_stokes_axisymmetric!(result, u_ops, qpinfo)
    u, ∇u, p = view(u_ops, 1:2), view(u_ops, 3:6), view(u_ops, 7)
    r = qpinfo.x[1]
    μ = qpinfo.params[1]
    # add Laplacian
    result[1] = μ / r * u[1] - p[1]
    result[2] = 0
    result[3] = μ * r * ∇u[1] - r * p[1]
    result[4] = μ * r * ∇u[2]
    result[5] = μ * r * ∇u[3]
    result[6] = μ * r * ∇u[4] - r * p[1]
    result[7] = -(r * (∇u[1] + ∇u[4]) + u[1])
    return nothing
end

function u!(result, qpinfo)
    x = qpinfo.x
    result[1] = x[1]
    result[2] = -2 * x[2]
    return nothing
end

function kernel_l2div(result, u_ops, qpinfo)
    u, divu = view(u_ops, 1:2), view(u_ops, 3)
    result[1] = (qpinfo.x[1] * divu[1] + u[1])^2
    return nothing
end


function main(; μ = 0.1, nrefs = 4, nonlinear = false, uniform = false, Plotter = nothing, kwargs...)

    # problem description
    PD = ProblemDescription()
    u = Unknown("u"; name = "velocity")
    p = Unknown("p"; name = "pressure")
    assign_unknown!(PD, u)
    assign_unknown!(PD, p)
    assign_operator!(PD, BilinearOperator(kernel_stokes_axisymmetric!, [id(u), grad(u), id(p)]; params = [μ], kwargs...)) #; jacobian = kernel_jacobian!))
    if nonlinear
        assign_operator!(PD, NonlinearOperator(kernel_convection!, [id(u)], [id(u), grad(u)]; bonus_quadorder = 1, kwargs...)) #; jacobian = kernel_jacobian!))
    end
    assign_operator!(PD, InterpolateBoundaryData(u, u!; regions = [3]))
    assign_operator!(PD, HomogeneousBoundaryData(u; regions = [4], mask = (1, 0, 1)))
    assign_operator!(PD, HomogeneousBoundaryData(u; regions = [1], mask = (0, 1, 1)))

    # grid
    if uniform
        xgrid = uniform_refine(grid_unitsquare(Triangle2D), nrefs)
    else
        xgrid = simplexgrid(
            Triangulate;
            points = [0 0 ; 5 0 ; 5 1 ; 0 1]',
            bfaces = [1 2 ; 2 3 ; 3 4 ; 4 1 ]',
            bfaceregions = [1, 2, 3, 4],
            regionpoints = [0.5 0.5;]',
            regionnumbers = [1],
            regionvolumes = [4.0^(-nrefs - 1)]
        )
    end

    # solve
    FES = [FESpace{H1P2B{2, 2}}(xgrid), FESpace{L2P1{1}}(xgrid)]
    sol = ExtendableFEM.solve(PD, FES; kwargs...)

    # compute divergence in cylindrical coordinates by volume integrals
    DivIntegrator = ItemIntegrator(kernel_l2div, [id(u), div(u)]; quadorder = 4, resultdim = 1)
    div_error = sqrt(sum(evaluate(DivIntegrator, sol)))
    @info "||div(u_h)|| = $div_error"

    # compute L2error
    function kernel_l2error(result, u_ops, qpinfo)
        u!(result, qpinfo)
        return result .= (result - u_ops) .^ 2
    end
    ErrorIntegratorExact = ItemIntegrator(kernel_l2error, [id(1)]; entities = ON_BFACES, regions = [3], quadorder = 4, kwargs...)
    error = evaluate(ErrorIntegratorExact, sol)
    L2error = sqrt(sum(view(error, 1, :)) + sum(view(error, 2, :)))
    @info "||u - u_h|| = $L2error"

    # plot
    plt = plot([id(u)], sol; Plotter = Plotter)

    return [div_error, L2error], plt
end

end # module

This page was generated using Literate.jl.