3D Incompressible Navier-Stokes: boundary condition problem

1

my error :

image
image926×67 1.92 KB

Can someone help me with this problem? Here are my full code file and link:
https://drive.google.com/drive/folders/1__zQLO3EFFaLeNT40qAm1GbA5wBKoBnf?usp=sharing
I think my problem is my boundary condition cannot be matched. Please try to help me solve this problem! I will appreciate it with you!

from dolfin import*
from VarMINT import *
from mshr import *
import numpy as np


####### Parameters #######

# Arguments are parsed from the command line, with some hard-coded defaults
# here.  See help messages for descriptions of parameters.

import argparse
parser = argparse.ArgumentParser()

parser.add_argument('--Re',dest='Re',default=10000.0,
                    help='Reynolds number.')
parser.add_argument('--k',dest='k',default=1,
                    help='Polynomial degree.')
parser.add_argument('--T',dest='T',default=3,
                    help='Length of time interval to consider.') 

args = parser.parse_args()

Re = Constant(float(args.Re))
k = int(args.k)
T = float(args.T)

####### Analysis #######

# Avoid overkill automatically-determined quadrature degree:
QUAD_DEG = 2*k
dx = dx(metadata={"quadrature_degree":QUAD_DEG})


# Domain:
import math

# mesh

mesh = Mesh()
filename = "cylinder.xdmf"
f = XDMFFile(filename)
f.read(mesh)


# mesh = RectangleMesh(Point(-math.pi,-math.pi),
                     # Point(math.pi,math.pi),
                     # Nel,Nel)

# Mixed velocity--pressure space:
V = equalOrderSpace(mesh,k=k)

# Solution and test functions:
up = Function(V)
u,p = split(up)
vq = TestFunction(V)
v,q = split(vq)

# Midpoint time integration:
N_STEPS = 1000 # Space--time quasi-uniformity
Dt = Constant(T/N_STEPS)
up_old = Function(V)
u_old,_ = split(up_old)
u_mid = 0.5*(u+u_old)
u_t = (u-u_old)/Dt

# Determine physical parameters from Reynolds number:
rho = Constant(1.0)
mu = Constant(1.0/Re)
nu = mu/rho

# Initial condition for the Taylor--Green vortex
x = SpatialCoordinate(mesh)

u_IC = as_vector((0*x[0], 0*x[1], 0*x[2]))

# Weak problem residual; note use of midpoint velocity:
F = interiorResidual(u_mid,p,v,q,rho,mu,mesh,
                     u_t=u_t,Dt=Dt,C_I=Constant(6.0*(k**4)),dx=dx)
                    
# Project the initial condition:
up_old.assign(project(as_vector((u_IC[0], u_IC[1], u_IC[2],Constant(0.0))),V))

# Same-velocity predictor:
up.assign(up_old)

inflow_profile = ('4.0*1.5*x[1]*(0.41 - x[1]) / pow(0.41, 2)', '0', '0')

#Define boundaries

inflow   = "near(x[0], 0)"
outflow  = "near(x[0], 12)"
walls    = "near(x[1], 0) || near(x[1], 5) || near(x[2], 0) || near(x[2], 3)"
cylinder = "on_boundary && x[0]>3 && x[0]<7 && x[1]>0.0 &&  x[1]<4.0 && x[2]>0.4 && x[2]<2.6"
                   
# Define boundary conditions
bcu_inflow = DirichletBC(V.sub(0), Expression(inflow_profile, degree=2), inflow)
bcu_walls = DirichletBC(V.sub(0), Constant((0, 0, 0)), walls)
bcu_cylinder = DirichletBC(V.sub(0), Constant((0, 0, 0)), cylinder)
bcp_outflow = DirichletBC(V.sub(1), Constant(0), outflow)
bcu = [bcu_inflow, bcu_walls, bcu_cylinder]
bcp = [bcp_outflow]
bcs = [bcu_inflow, bcu_walls, bcu_cylinder, bcp_outflow]


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

# Time stepping loop:
for step in range(0,N_STEPS):
    print("======= Time step "+str(step+1)+"/"+str(N_STEPS)+" =======")
    solve(F==0,up,bcs=bcs)
    up_old.assign(up)
    
    # Extract solution components
    u, p = up.split() # seperate functions
    u.rename("u", "velocity")
    p.rename("p", "pressure")
    
    # Update previous solution
    if (step % (N_STEPS/100) == 0):
        ufile << u
        # pfile << p
    	# Save to file
2

Please do not post the same questions twice.
Also, as please reduce the problem to a minimal example. As it is the assignment of the boundary conditions that is not working, you can remove anything not related to the boundary conditions.

As your mesh is based on a GMSH file, you should use Physical Markers in gmsh to assign your boundary conditions. Then you will not have issues with them not being found through your C++ conditionals.

4

Instead of assigning a single physical tag to the surfaces (19 in your case), create a physical surface for each of the groups that corresponds to a boundary condition. Then you can load these into dolfin using meshio, as explained at: Mesh generation and conversion with GMSH and PYGMSH | Jørgen S. Dokken

6

The point is that you should write a triangular mesh (as shown in the tutorial) alongside your tetrahedral mesh.
You then read in the triangular mesh and load it into a MeshValueCollection, as described here:

Then, you can use DirichletBC in the way shown here: