245 : Flow around a cylinder

(source code)

This example solves the DFG Navier-Stokes benchmark problem

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

on a rectangular 2D domain with a circular obstacle, see here for details.

This script demonstrates the employment of external grid generators and the computation of drag and lift coefficients.

Note: This example needs the additional packages Triangulate and SimplexGridFactory to generate the mesh.

module Example245_NSEFlowAroundCylinder

using ExtendableFEM
using Triangulate
using SimplexGridFactory
using ExtendableGrids
using GridVisualize
using LinearAlgebra

# inlet data for Karman vortex street example
# as in DFG benchmark 2D-1 (Re = 20, laminar)
const umax = 0.3
const umean = 2 // 3 * umax
const L, W, H = 0.1, 2.2, 0.41
function inflow!(result, qpinfo)
    x = qpinfo.x
    result[1] = 4 * umax * x[2] * (H - x[2]) / (H * H)
    result[2] = 0.0
    return nothing
end

# hand constructed identity matrix for kernel to avoid allocations
const II = [1 0; 0 1]


# Example of a kernel using tensor_view() function to allow for an operator
# based style of writing the semilinear form.
# For comparison we also provide the kernel_nonlinear_flat! function below
# that uses a component-wise style of writing the semilinear form.

#
# the scalar product ``(\nabla v, \mu \nabla u - pI)`` will be evaluated
# so in general `a = b` corresponds to ``(a,b)``.

# Note that the order of vector entries between the kernel and the call to
# NonlinearOperator have to match.
function kernel_nonlinear!(result, u_ops, qpinfo)
    # Shape values of vectorial u are starting at index 1
    # view as 1-tensor(vector) of length dim=2 in 2D
    u = tensor_view(u_ops, 1, TDVector(2))
    v = tensor_view(result, 1, TDVector(2))
    # gradients of vectorial u are starting at index 3
    # view as 2-tensor of size 2x2 in 2D
    ∇u = tensor_view(u_ops, 3, TDMatrix(2))
    ∇v = tensor_view(result, 3, TDMatrix(2))
    # values of scalar p are starting at index 7
    # view as 0-tensor (single value)
    p = tensor_view(u_ops, 7, TDScalar())
    q = tensor_view(result, 7, TDScalar())
    # get viscosity at current quadrature point
    μ = qpinfo.params[1]
    # Note that all operators should be element-wise to avoid allocations
    # `(v,u⋅∇u) = (v,∇u^T⋅u)`
    tmul!(v, ∇u, u)
    # `(∇v,μ∇u-pI)`
    ∇v .= μ .* ∇u .- p[1] .* II
    # `(q,-∇⋅u)`
    q[1] = -dot(∇u, II)
    return nothing
end


# everything is wrapped in a main function
function main(; Plotter = nothing, μ = 1.0e-3, maxvol = 1.0e-3, reconstruct = true, parallel = false, npart = 8, kwargs...)

    # load grid (see function below)
    xgrid = make_grid(W, H; n = Int(ceil(sqrt(1 / maxvol))), maxvol = maxvol)
    if parallel
        xgrid = partition(xgrid, RecursiveMetisPartitioning(npart = npart))
    end

    # problem description
    PD = ProblemDescription()
    u = Unknown("u"; name = "velocity")
    p = Unknown("p"; name = "pressure")
    id_u = reconstruct ? apply(u, Reconstruct{HDIVRT1{2}, Identity}) : id(u)

    assign_unknown!(PD, u)
    assign_unknown!(PD, p)
    assign_operator!(PD, NonlinearOperator(kernel_nonlinear!, [id_u, grad(u), id(p)]; params = [μ], parallel = parallel, kwargs...))
    assign_operator!(PD, InterpolateBoundaryData(u, inflow!; regions = 4))
    assign_operator!(PD, HomogeneousBoundaryData(u; regions = [1, 3, 5]))

    # P2-bubble + reconstruction operator
    FETypes = [H1P2B{2, 2}, H1P1{1}]

    # generate FESpaces and Solution vector
    FES = [FESpace{FETypes[1]}(xgrid), FESpace{FETypes[2]}(xgrid; broken = true)]

    # solve
    sol = solve(PD, FES; maxiterations = 50, target_residual = 1.0e-10)

    # postprocess : compute drag/lift (see function below)
    draglift = get_draglift(sol, μ; parallel = parallel, kwargs...)
    pdiff = get_pressure_difference(sol)
    println("[drag, lift] = $draglift")
    println("p difference = $pdiff")

    # plots via GridVisualize
    plt = GridVisualizer(; Plotter = Plotter, layout = (4, 1), clear = true, size = (800, 1200))
    gridplot!(plt[1, 1], xgrid, cellcoloring = :partitions, linewidth = 1)
    gridplot!(plt[2, 1], xgrid, cellcoloring = :partitions, linewidth = 1, xlimits = [0, 0.3], ylimits = [0.1, 0.3])
    scalarplot!(plt[3, 1], xgrid, nodevalues(sol[u]; abs = true)[1, :])
    vectorplot!(plt[3, 1], xgrid, eval_func_bary(PointEvaluator([id(u)], sol)), rasterpoints = 20, clear = false)
    scalarplot!(plt[4, 1], xgrid, view(nodevalues(sol[p]), 1, :), levels = 11, title = "p_h")

    return [draglift[1], draglift[2], pdiff[1]], plt
end

function get_pressure_difference(sol::FEVector)
    xgrid = sol[2].FES.xgrid
    PE = PointEvaluator([id(2)], sol)
    p_left = zeros(Float64, 1)
    x1 = [0.15, 0.2]
    p_right = zeros(Float64, 1)
    x2 = [0.25, 0.2]
    evaluate!(p_left, PE, x1)
    evaluate!(p_right, PE, x2)
    @show p_left, p_right
    return p_left - p_right
end

function get_draglift(sol::FEVector, μ; parallel = false, kwargs...)

    # this function is interpolated for drag/lift test function creation
    function DL_testfunction(component)
        return function closure(result, qpinfo)
            x = qpinfo.x
            fill!(result, 0)
            return if sqrt((x[1] - 0.2)^2 + (x[2] - 0.2)^2) <= 0.06
                result[component] = 1
            end
        end
    end

    # drag lift calculation by testfunctions
    function draglift_kernel(result, input, qpinfo)
        # input = [ u, grad(u), p , v , grad(v)]
        #         [1:2,  3:6,   7 ,8:9,  10:13 ]
        result[1] = μ * (input[3] * input[10] + input[4] * input[11] + input[5] * input[12] + input[6] * input[13])
        result[1] += (input[1] * input[3] + input[2] * input[4]) * input[8]
        result[1] += (input[1] * input[5] + input[2] * input[6]) * input[9]
        result[1] -= input[7] * (input[10] + input[13])
        result[1] *= -(2 / (umean^2 * L))
        return nothing
    end
    DLIntegrator = ItemIntegrator(draglift_kernel, [id(1), grad(1), id(2), id(3), grad(3)]; quadorder = 4, parallel = parallel, kwargs...)

    # test for drag
    TestFunction = FEVector(sol[1].FES; name = "drag/lift testfunction")
    interpolate!(TestFunction[1], ON_BFACES, DL_testfunction(1))
    drag = sum(evaluate(DLIntegrator, [sol[1], sol[2], TestFunction[1]]))

    # test for lift
    interpolate!(TestFunction[1], ON_BFACES, DL_testfunction(2))
    lift = sum(evaluate(DLIntegrator, [sol[1], sol[2], TestFunction[1]]))

    return [drag, lift]
end

# grid generator script using SimplexGridBuilder/Triangulate
function make_grid(W, H; n = 20, maxvol = 0.1)
    builder = SimplexGridBuilder(Generator = Triangulate)
    function circlehole!(builder, center, radius; n = 20)
        points = [point!(builder, center[1] + radius * sin(t), center[2] + radius * cos(t)) for t in range(0, 2π, length = n)]
        for i in 1:(n - 1)
            facet!(builder, points[i], points[i + 1])
        end
        facet!(builder, points[end], points[1])
        return holepoint!(builder, center)
    end
    p1 = point!(builder, 0, 0)
    p2 = point!(builder, W, 0)
    p3 = point!(builder, W, H)
    p4 = point!(builder, 0, H)

    # heuristic refinement around cylinder
    refine_radius = 0.25
    maxrefinefactor = 1 // 20
    function unsuitable(x1, y1, x2, y2, x3, y3, area)
        if area > maxvol * min(max(4 * maxrefinefactor, abs((x1 + x2 + x3) / 3 - 0.2)), 1 / maxrefinefactor)
            return true
        end
        dist = sqrt(((x1 + x2 + x3) / 3 - 0.2)^2 + ((y1 + y2 + y3) / 3 - 0.2)^2) - 0.05
        myarea = dist < refine_radius ? maxvol * max(maxrefinefactor, 1 - (refine_radius - dist) / refine_radius) : maxvol
        if area > myarea
            return true
        else
            return false
        end
    end

    facetregion!(builder, 1)
    facet!(builder, p1, p2)
    facetregion!(builder, 2)
    facet!(builder, p2, p3)
    facetregion!(builder, 3)
    facet!(builder, p3, p4)
    facetregion!(builder, 4)
    facet!(builder, p4, p1)
    facetregion!(builder, 5)
    circlehole!(builder, (0.2, 0.2), 0.05, n = n)

    return simplexgrid(builder, maxvolume = 16 * maxvol, unsuitable = unsuitable)
end

end

This page was generated using Literate.jl.