Integrating over an interior surface

variational formulation
9

Consider the following problem for a circular object:

from fenics import *
from mshr import *
import numpy as np

def compute_circulation(N):
    # Create mesh
    length = 4.0
    diameter = 1.0
    p0 = Point(np.array([0.0, 0.0]))
    p1 = Point(np.array([length, diameter]))
    r = 0.2
    obj = Circle(Point(length/2, diameter/2), r)
    domain = Rectangle(p0, p1) - obj
    mesh = generate_mesh(domain, N)
    V = FunctionSpace(mesh, 'P', 1)
    # Define boundary condition
    u_D = Expression('230*x[1]', degree=1)

    # Making a mark dictionary
    mark = {"generic": 0,"lower_wall": 1,"upper_wall": 2,"left": 3,"right": 4, "airfoil": 5 }
    subdomains = MeshFunction("size_t", mesh, 1)
    subdomains.set_all(mark["generic"])
    class Left(SubDomain):
        def inside(self, x, on_boundary):
            return on_boundary and near(x[0], 0)
    class Right(SubDomain):
        def inside(self, x, on_boundary):
            return on_boundary and near(x[0], length)
    class UpperWall(SubDomain):
        def inside(self, x, on_boundary):
            return on_boundary and near(x[1], diameter)
    class LowerWall(SubDomain):
        def inside(self, x, on_boundary):
            return on_boundary and near(x[1], 0)
    class Airfoil(SubDomain):
        def inside(self, x, on_boundary):
            return on_boundary and (x[0]-length/2)**2+(x[1]-diameter/2)**2<1.1*r**2
    File("mesh.pvd") << mesh
    # Marking the subdomains
    left = Left()
    left.mark(subdomains, mark["left"])
    right = Right()
    right.mark(subdomains, mark["right"])
    upper_wall = UpperWall()
    upper_wall.mark(subdomains, mark["upper_wall"])
    lower_wall = LowerWall()
    lower_wall.mark(subdomains, mark["lower_wall"])
    airfoil = Airfoil()
    airfoil.mark(subdomains, mark["airfoil"])

    bc_lower_wall = DirichletBC(V, 0, subdomains, mark["lower_wall"])
    bc_upper_wall = DirichletBC(V, u_D((0,diameter)), subdomains, mark["upper_wall"])
    bc_left = DirichletBC(V, u_D, subdomains, mark["left"])
    bc_right = DirichletBC(V, u_D, subdomains, mark["right"])
    bc_airfoil = DirichletBC(V, u_D((2,diameter/2)), subdomains, mark["airfoil"])
    bcs = [bc_lower_wall, bc_upper_wall, bc_left, bc_right, bc_airfoil]
    # Define variational problem
    u = TrialFunction(V)
    v = TestFunction(V)
    f = Constant(0)
    a = dot(grad(u), grad(v))*dx
    L = f*v*dx
    # Compute solution
    u = Function(V)
    solve(a == L, u, bcs)
    W = VectorFunctionSpace(mesh, 'P', 1)
    vx = project(u.dx(1), V)
    vy = project(-u.dx(0), V)
    vvec = project(as_vector([vx, vy]),W)

    ds = Measure("ds", domain=mesh, subdomain_data=subdomains)
    GammaP = ds(5)
    n = FacetNormal(mesh)
    tangent = as_vector([n[1], -n[0]])
    L1 = (dot(vvec, tangent))*GammaP

    # Project tangent to CG 1 space for visualization
    VV = VectorFunctionSpace(mesh, "CG", 1)
    q,r = TrialFunction(VV), TestFunction(VV)
    aV = inner(q, r)*ds(5)
    lV = inner(tangent, r)*ds(5)
    tang_proj = Function(VV)
    AV = assemble(aV, keep_diagonal=True)
    AV.ident_zeros()
    LV = assemble(lV)
    solve(AV, tang_proj.vector(), LV)
    File("tang.pvd") << tang_proj
    File("circ.pvd") << vvec

    circulation = assemble(L1)
    print(N, circulation, assemble(inner(vvec, as_vector((1,0)))*ds(5)))


compute_circulation(50)
compute_circulation(100)
compute_circulation(200)
compute_circulation(400)

yielding:

fenics@3a01410489e7:~/shared$ python3 circulation.py 
Solving linear variational problem.
50 -2.814902850116198 300.78413974646304
100 -0.35969509592035376 314.1475947113785
200 -0.3435394737202848 323.7011340127153
400 0.03041307517778638 328.270430678979

Here I have made the circulation a function of the mesh resolution (for a circular object). As you can observe, on a coarse mesh, the mesh circulation is far from 0 on the coarsest meshes. Also note that last value is simply the integral of the x-component of the velocity field, is of quite a large magnitude (and the flow field around the obstacle does not satisfy a no penetration constraint.

Therefore I would suggest using a finer mesh.

Note that I have also added a projection of the tangent vector to a CG 1 space, such that it can be visualized with either matplotlib or Paraview.

The accuracy of the circulation integral should also be greatly improved by using a second order mesh. There is limited support for higher order geometries in dolfin, while this has been one of the main focuses of dolfinx

2 Likes