Obtain velocity normal to boundary

1

Hi everyone. I hope you can guide me to further my studies. I’m working on a 2D free surface problem with a sine wave as initial shape for the top boundary. I solve the 2D Navier-Stokes equations using a mixed formulation along with 2 boundary conditions for the free surface (namely a stress condition “n.T = sigma*kappa” and a kinematic condition for the velocity at the interface “u.n = V.n”.

I then take the solution of the system (i.e. the velocity vector) and project it onthe surface coordinate system (normal n, tangential t). I calculate the surface displacement with the velocity normal to the surface, in this manner:

u_normal = inner(n, u)*n

However the type of u_normal is a ComponentTensor, and I can’t figure out how to access the indices of this object. Is there a more intelligent way to perform such operation? I attach below a MWE.

from dolfin import *
import numpy as np
import matplotlib.pyplot as plt

class Top(SubDomain):
    def __init__(self):
        SubDomain.__init__(self)
    def inside(self, x, on_boundary):
        return near(x[1], 1)

# Arbitrary time step to get displacement
dt = 0.01
# Create mesh and boundary mesh
mesh =  UnitSquareMesh(20,20)
boundary_mesh = BoundaryMesh(mesh, "exterior", True) 

# Mesh markers    mesh_markers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)    boundary_markers = MeshFunction("size_t", boundary_mesh, 0)

# Mark mesh
mesh_markers.set_all(0)
Top().mark(mesh_markers, 1)

# mark boundary_mesh
boundary_markers.set_all(0)
Top().mark(boundary_markers, 1)

# Create function space for bulk dynamics
V = VectorFunctionSpace(mesh, "CG", 2)
u = Function(V)

# Assign random values to bulk velocity    u.vector().set_local(np.random.rand(u.vector().local_size()))
# Create function space for interface dynamics
D= VectorFunctionSpace(boundary_mesh, "CG", 1)
displacement = Function(D)
displacement.set_allow_extrapolation(True)

# Get normal
n = FacetNormal(mesh)

# Obtain velocity normal vector 
#! NOT WORKING
u_normal = n * inner(n, u)

# Interpolate mesh values into boundary function
u_boundary = Function(D)
u_boundary.set_allow_extrapolation(True)

# Interpolate velocity normal vector into boundary function
u_boundary.interpolate(u_normal)
# Kinematic boundary condition for free surface leads to boundary displacement 
# displacement = u_boundary*dt    top_disp = DirichletBC(D, dt*u_boundary, boundary_markers, 1)

# apply DirichletBC into displacement function
top_disp.apply(displacement.vector())

# move mesh
ALE.move(boundary_mesh, displacement)
ALE.move(mesh, boundary_mesh)
mesh.bounding_box_tree().build(mesh)

plot(mesh); plt.show()
2

As the FacetNormal as implied by its name, is defined on each Facet, and not for each degree of freedom, you have to approximate the normal at the degrees of freedom (for instance at corners).

For your particular case, the problem is that you are trying to interpolate an ComponentTensor, that is an expression using ufl syntax. To work around this, I would suggest using an UserExpression as follows:

n = FacetNormal(mesh)
# Interpolate mesh values into boundary function
u_boundary = Function(D)
u_boundary.set_allow_extrapolation(True)

# Approximate facet normal in a suitable space using projection
n = FacetNormal(mesh)
V = VectorFunctionSpace(mesh, "CG", 2)
u_ = TrialFunction(V)
v_ = TestFunction(V)
a = inner(u_,v_)*ds
l = inner(n, v_)*ds
A = assemble(a, keep_diagonal=True)
L = assemble(l)

A.ident_zeros()
nh = Function(V)

solve(A, nh.vector(), L)

class normal_u(UserExpression):
    def __init__(self, **kwargs):
        super().__init__(**kwargs)
        
    def eval(self, values, x):
        n_eval = nh(x)
        u_eval = u(x)
        un = (u_eval[0]*n_eval[0] + u_eval[1]*n_eval[1])
        values[0] = un * n_eval[0]
        values[1] = un * n_eval[1]

    def value_shape(self):
        return (2,)

u_boundary.interpolate(normal_u())
3

Thank you for the reply! Your solution solved my problem.

4

Hi @heitorvc , If you have computed the shear stresses for the 2D Navier-Stokes solution, can you please share the code