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:
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?