T = 5
num_steps = 5000
dt = T / num_steps
mu = 0.01
rho = 1
channel = Rectangle(Point(0, 0), Point(2.2, 0.41))
# gen_mesh is a self-written function that returns a Polygon object
cylinder = gen_mesh(c=np.array([1.5, 0.05]), n=5, scale=0.3, rad=0.2, edgy=0.05)
domain = channel - cylinder
mesh = generate_mesh(domain, 64)
# Define function spaces
V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 1)
# Define boundaries
inflow = 'near(x[0], 0)'
outflow = 'near(x[0], 2.2)'
walls = 'near(x[1], 0) || near(x[1], 0.41)'
cylinder = 'on_boundary && x[0] > 1 && x[0] < 2 && x[1] > 0.1 && x[1] < 0.4'
# Define inflow profile
inflow_profile = ('4.0 * 1.5 * x[1] * (0.41 - x[1]) / pow(0.41, 2)', '0')
# Define boundary conditions
bcu_inflow = DirichletBC(V, Expression(inflow_profile, degree=2), inflow)
bcu_walls = DirichletBC(V, Constant((0, 0)), walls)
bcu_cylinder = DirichletBC(V, Constant((0, 0)), cylinder)
bcp_outflow = DirichletBC(Q, Constant(0), outflow)
bcu = [bcu_inflow, bcu_walls, bcu_cylinder]
bcp = [bcp_outflow]
I am looking to solve a channel flow problem with an obstacle placed in the channel. The boundary conditions are defined as above based on the example code given in the tutorial. However, the boundary conditions are violated as shown in the velocity distribution plot below after one time step, is there an error with my code? I suspect the problem comes from
cylinder = 'on_boundary && x[0] > 1 && x[0] < 2 && x[1] > 0.1 && x[1] < 0.4'
but I am not entirely sure. Any help is appreciated, thanks.
