Normal velocity (vector) Dirichlet boundary condition used for a surface of 3d mesh

[xirancao1]

When I used a normal velocity (vector) Dirichlet boundary condition for a surface of 3d mesh as inlet condition, the velocity on edges near the wall was not a constant? The plot is as follows. Can anyone help me ?

And my code of the mesh generation is below, the .py file was uploded to Github as pipe.py, the .stl file used in the script is pipe.stl, the generated .msh file is pipe.msh:

GitHub - caoxiran/fenicsxx

import gmsh
import sys
import os
import math
import meshio
gmsh.initialize(sys.argv)

# merge STL, create surface patches that are reparametrizable (so we can remesh
# them) and compute the parametrizations
path = os.path.dirname(os.path.abspath(__file__))
gmsh.merge(os.path.join(path, 'pipe.stl'))
gmsh.model.mesh.classifySurfaces(math.pi, True, True)

#classifysurface:Classify ("color") the current surface mesh based on an angle threshold (the first argument, in radians), and create new discrete surfaces, curves and points accordingly. If the second argument is set, also create discrete curves on the boundary if the surface is open. If the third argument is set, create edges and surfaces than can be reparametrized with CreateGeometry. The last optional argument sets an angle threshold to force splitting of the generated curves
gmsh.model.mesh.createGeometry()
# Create a geometry for discrete entities (represented solely by a mesh, without an underlying CAD description) in the current model, i.e. create a parametrization for discrete curves and surfaces, assuming that each can be parametrized with a single map. If no entities are given, create a geometry for all discrete entities. 
#extrude = False
#extrude = False


#Get all the entities in the current model. A model entity is represented by two integers: its dimension (dim == 0, 1, 2 or 3) and its tag (its unique, strictly positive identifier). If dim is >= 0, return only the entities of the specified dimension (e.g. points if dim == 0). The entities are returned as a vector of (dim, tag) pairs. 
    # get "top" surfaces created by extrusion
gmsh.model.mesh.optimize('Relocate2D',True)
sur=gmsh.model.getEntities(2)
surface = [s[1] for s in sur]
top_ent = gmsh.model.getEntities(2)
top_surf = [s[1] for s in top_ent]
print(top_surf)
# get boundary of top surfaces, i.e. boundaries of holes
gmsh.model.geo.synchronize()
bnd_ent = gmsh.model.getBoundary(top_ent)
bnd_curv = [c[1] for c in bnd_ent]
print(top_surf)
print(bnd_curv)
# create plane surfaces filling the holes
loops = gmsh.model.geo.addCurveLoops(bnd_curv)
print(top_surf)
print(bnd_curv)
print(loops)
gmsh.model.mesh.optimize('Relocate2D',True)
for l in loops:
    top_surf.append(gmsh.model.geo.addPlaneSurface([l]))
sur_boundary=gmsh.model.getEntities(2)
surface_boundary = [s[1] for s in sur_boundary]
print(top_surf)
print(bnd_curv)
print(loops)
print(top_surf)
print(surface_boundary)
gmsh.model.mesh.optimize('Relocate2D',True)

in_marker=1
out_marker=2
wall_marker=3
domain_marker=4
gmsh.model.geo.addPhysicalGroup(2,[4],in_marker)
gmsh.model.setPhysicalName(2,in_marker,'in')
gmsh.model.geo.addPhysicalGroup(2,[5],out_marker)
gmsh.model.setPhysicalName(2,out_marker,'out')
gmsh.model.geo.addPhysicalGroup(2,surface,wall_marker)
gmsh.model.setPhysicalName(2,wall_marker,'wall')
#gmsh.model.geo.addPhysicalGroup(2,[177],6)
#gmsh.model.setPhysicalName(2,6,'boundary')
gmsh.model.geo.addPhysicalGroup(3,[1],domain_marker)
gmsh.model.setPhysicalName(3,domain_marker,'domain')
    # create the inner volume
gmsh.model.geo.addVolume([gmsh.model.geo.addSurfaceLoop(top_surf)])
gmsh.model.geo.synchronize()

# use MeshAdapt for the resulting not-so-smooth parametrizations
gmsh.option.setNumber('Mesh.Algorithm', 1)
gmsh.option.setNumber('Mesh.MeshSizeFactor', 0.3)
gmsh.model.mesh.generate()
gmsh.write('pipe.msh')
#Launch the GUI to see the results:
if '-nopopup' not in sys.argv:
    gmsh.fltk.run()

gmsh.finalize()
msh=meshio.read("tipe.msh")
meshio.write("pipe.xdmf", meshio.Mesh(points=msh.points, cells={"tetra": msh.cells["tetra"]}))
#boundary face
meshio.write("mf_pipe.xdmf", meshio.Mesh(points=msh.points, cells={"triangle": msh.cells["triangle"]},cell_data={"triangle": {"name_to_read": msh.cell_data["triangle"]["gmsh:physical"]}}))
#boundary cell
meshio.write("cf_pipe.xdmf", meshio.Mesh(points=msh.points, cells={"tetra": msh.cells["tetra"]},cell_data={"tetra": {"name_to_read":msh.cell_data["tetra"]["gmsh:physical"]}}))

The Navier-Stokes equations was solved in demo_navier-stokes_3dtest.py, the code is as follows:

"""This demo program solves the incompressible Navier-Stokes equations
on an L-shaped domain using Chorin's splitting method."""

# Copyright (C) 2010-2011 Anders Logg
#
# This file is part of DOLFIN.
#
# DOLFIN is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# DOLFIN is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
#
# Modified by Mikael Mortensen 2011
#
# First added:  2010-08-30
# Last changed: 2011-06-30

# Begin demo


import matplotlib.pyplot as plt
from dolfin import *

# Print log messages only from the root process in parallel
parameters["std_out_all_processes"] = False;

# Load mesh from file
#mesh = Mesh("lshape.xml.gz")
mesh = Mesh()
with XDMFFile("pipe.xdmf") as infile:
    infile.read(mesh)
##
mvc = MeshValueCollection("size_t", mesh, 2) 
with XDMFFile("mf_pipe.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
boundary_markers=MeshFunction("size_t",mesh,mvc)
##
mvc = MeshValueCollection("size_t", mesh, 3)
with XDMFFile("cf_pipe.xdmf") as infile:
    infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
# Define function spaces (P2-P1)
V = VectorFunctionSpace(mesh, "Lagrange", 3)
Q = FunctionSpace(mesh, "Lagrange", 1)

# Define trial and test functions
u = TrialFunction(V)
p = TrialFunction(Q)
v = TestFunction(V)
q = TestFunction(Q)
in_marker=2
out_marker=1
wall_marker=3
domain_marker=4
# Set parameter values
dt = 0.001
T = 1
#nu = 0.01
nu=0.0000033018868
# Define time-dependent pressure boundary condition
#p_in = Expression("sin(3.0*t)", t=0.0, degree=2)

# Define boundary conditions
#noslip  = DirichletBC(V, (0, 0),           "on_boundary && \                       (x[0] < DOLFIN_EPS | x[1] < DOLFIN_EPS | \                       (x[0] > 0.5 - DOLFIN_EPS && x[1] > 0.5 - DOLFIN_EPS))")
#inflow  = DirichletBC(Q, p_in, "x[1] > 1.0 - DOLFIN_EPS")
#outflow = DirichletBC(Q, 0, "x[0] > 1.0 - DOLFIN_EPS")
#bcu = [noslip]
#bcp = [inflow, outflow]
bcu_in=DirichletBC(V,Constant((0.182866,-0.03869,-0.071156)),boundary_markers,in_marker)

bcp_out=DirichletBC(Q,Constant(666.61),boundary_markers,out_marker)
bcu_wall=DirichletBC(V,Constant((0.0,0.0,0.0)),boundary_markers,wall_marker)
bcu=[bcu_in,bcu_wall]
bcp=[bcp_out]
# Create functions
u0 = Function(V)
u1 = Function(V)
p1 = Function(Q)

# Define coefficients
k = Constant(dt)
f = Constant((0, 0,0))

# Tentative velocity step
F1 = (1/k)*inner(u - u0, v)*dx + inner(grad(u0)*u0, v)*dx + \
     nu*inner(grad(u), grad(v))*dx - inner(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)

# Pressure update
a2 = inner(grad(p), grad(q))*dx
L2 = -(1/k)*div(u1)*q*dx

# Velocity update
a3 = inner(u, v)*dx
L3 = inner(u1, v)*dx - k*inner(grad(p1), v)*dx

# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

# Use amg preconditioner if available
prec = "amg" if has_krylov_solver_preconditioner("amg") else "default"

# Use nonzero guesses - essential for CG with non-symmetric BC
parameters['krylov_solver']['nonzero_initial_guess'] = True

# Create files for storing solution
ufile = File("results/velocity3d.pvd")
pfile = File("results/pressure3d.pvd")

# Time-stepping
t = dt
while t < T + DOLFIN_EPS:
    print(t)
    # Update pressure boundary condition
    #p_in.t = t
    begin("Computing tentative velocity")
    # Compute tentative velocity step
    b1 = assemble(L1)
    [bc.apply(A1, b1) for bc in bcu]
    solve(A1, u1.vector(), b1, "bicgstab", "default")
    end()
    begin("Computing pressure correction")
    # Pressure correction
    b2 = assemble(L2)
    [bc.apply(A2, b2) for bc in bcp]
    [bc.apply(p1.vector()) for bc in bcp]
    solve(A2, p1.vector(), b2, "bicgstab", prec)
    end()
    begin("Computing velocity correction")
    # Velocity correction
    b3 = assemble(L3)
    [bc.apply(A3, b3) for bc in bcu]
    solve(A3, u1.vector(), b3, "bicgstab", "default")
    end()
    # Save to file
    ufile << u1
    pfile << p1

    # Move to next time step
    u0.assign(u1)
    t += dt

I really know why the velocity boundary condition was different at the edges of the surface, Could anyone help me solve the issue?

[dokken]

I do not understand your issue.
are you confused why the solution gradually goes to zero at the edge of the circular inlet?
This is due to the fact that you set a zero bc on the cylindrical wall, which is set after you set the inlet condition, efficiently making it zero at those dofs.

[xirancao1]

Thank you for your reply! If the order of assigning boundary conditions between the wall and the inlet is changed, there will be no gradual decrease in the circular inlet surface？ I will try to verify the exchange order.

I also have another question, it is that does the grid drawn by Gmsh have no units when the. msh file is opened? If the coordinate units of grid points do not match the input parameter units in the equation, do we need to change the weak solution form of the equation while changing the boundary conditions? For example, if the coordinates in the NS equation are changed to mm, do the corresponding velocity and pressure boundary conditions need to be converted, and do the weak solution forms need to be changed?

[dokken]

@xirancao1

Thank you for your reply! If the order of assigning boundary conditions between the wall and the inlet is changed, there will be no gradual decrease in the circular inlet surface？ I will try to verify the exchange order.

If you change the order of application of your bcs, you should see no gradual decrease at the inlet, but you would see an instantaneous narrowing of the flow after the first element.

@xirancao1

I also have another question, it is that does the grid drawn by Gmsh have no units when the. msh file is opened? If the coordinate units of grid points do not match the input parameter units in the equation, do we need to change the weak solution form of the equation while changing the boundary conditions? For example, if the coordinates in the NS equation are changed to mm, do the corresponding velocity and pressure boundary conditions need to be converted, and do the weak solution forms need to be changed?

DOLFIN/DOLFINx does not assume any specific units. If you read in a mesh defined in millimeter, your should ensure that the units of your problem makes sense.

[xirancao1]

Thank you very much!