Neumann BC with point integral fundamentals in fenics

variational formulation
1

Hi everyone,

I’m trying to impose a Laplace-Young stress jump condition on a moving boundary (free surface) from a 2D domain. The expression is integrated by parts, resulting in

image

I’m struggling to evaluate the term v\cdot t \Big |^{S_i}_{S_f}. I obtain the tangent vector by rotating the normal vector as

import ufl
n = FacetNormal(mesh)
t_hat = ufl.perp(n)

I also interpolate the tangent into the mesh dofs. Currently I mark the moving boundary with flag “surface” and then mark its first and last facets with a different flag:

class Surf(SubDomain): 
    def __init__(self):
        SubDomain.__init__(self)

    def inside(self, x, on_boundary):
        return near(x[1], H)

class SL(SubDomain):
    def __init__(self):
        SubDomain.__init__(self)

    def inside(self, x, on_boundary):
        return near(x[1], H) and between(x[0], (0, xsl))

class SurfR(SubDomain):
    def __init__(self):
        SubDomain.__init__(self)

    def inside(self, x, on_boundary):
        return near(x[1], H) and between(x[0], (xsr, L))

in which xsl and xsr are manually defined (not the most intelligent way) in order to mark the respective facet of the boundary. I don’t know if this the proper way to do it. The red facet in figure below represents the surface left extremity.

image

I assure that tangent is nonzero on the surface extremities. Then, I add the stress into the variational problem by the following term (Neumann BC):

Uel = VectorElement('Lagrange', mesh.ufl_cell(), 2)
Pel = FiniteElement('Lagrange', mesh.ufl_cell(), 1)
UPel = MixedElement([Uel,Pel])
W = FunctionSpace(mesh, UPel)
(v, q) = TestFunctions(W)

I_S = Identity(mesh.geometry().dim()) - outer(n,n)

gradsv = I_S*grad(v)*t_hat

b_int = sigma*(- inner(gradsv, t_hat)*ds(surface_marker) +\
        inner(t_hat,v)*ds(surf_right) +\
      - inner(t_hat,v)*ds(surf_left) )

However I am obtaining a solution that blows at the surface extremities. I’m working on a slender domain, so it is very difficult for me to provide visual description. Figure above represents the pressure field obtained by the aforementioned approach.

image
image1280×320 26.6 KB

The only way to visualize the result is through a log scale, so the regions with pressure equal unity are negative pressure. The solution seems to make physical sense inside the domain, but the simulation does not converge after the first iteration due to the blowup of the left and right extremities. The expected result should be something like

image
image1233×105 16.9 KB

My question is: is there something fundamentally wrong with my approach? Does this occur because I’m using a test function from a mixed function space? Do I need to invest in learning dolfinx in order to work with calculations from different topologies together?

From the FEM perspective, the point integral terms should be added into the residuals of the respective elements, but I still didnt figure out a way to do it. I am working to simplify the problem so I can provide a MWE. However, the code is found here. Many thanks.