Unphysical result in using Quad grid in Stokes problem with pressure Direchlet boundary condition at inlet&outlet


131
views
0
3 months ago by
Dear FEniCS user,

I'm now trying to use Quad grid in a channel flow problem with Stokes equation, where pressure Direchlet boundary condition specified at inlet&outlet. The same problem I solved with triangle mesh is fine.
But now the result show that the pressure varying from p_in to p_out only at a very small region (1/10 of the streamwise) near inlet, and the rest part p=p_out.

But using Quad grid with velocity Direchlet boundary conditon does NOT have this problem.

My dolfin-version is: 2017.2.0.

The code is listed below, please let me know if you also have this trouble or you know the possible reason. Any comments or discussions are greatly appreciated!
from dolfin import *

mesh = UnitSquareMesh.create(30, 30, CellType.Type_quadrilateral)

mark = {"generic": 0,
"wall": 1,
"left": 2,
"right": 3 }

subdomains = MeshFunction("size_t", mesh, 1)
subdomains.set_all(mark["generic"])

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0)

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 1)

#Not sure how to define a class for both up and down boundary
class Wall(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 0) or near(x[1], 1)

left = Left()
left.mark(subdomains, mark["left"])

right = Right()
right.mark(subdomains, mark["right"])

wall = Wall()
wall.mark(subdomains, mark["wall"])

# Define function spaces
V = VectorElement("Lagrange", mesh.ufl_cell(), 2)
Q = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, V*Q)  # T-H element

# Define variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)

# Volume integration
dx = Measure("dx", domain=mesh, subdomain_data=subdomains)
# Surface integration
ds = Measure("ds", domain=mesh, subdomain_data=subdomains)
# Surface normal
n = FacetNormal(mesh)

# Pressures. First define the numbers (for later use):
p_left = -10000.0
p_right = 0.0
# ...and then the DOLFIN constants:
pressure_left = Constant(p_left)
pressure_right = Constant(p_right)
# Body force:
force = Constant((0.0, 0.0))

a = inner(grad(u), grad(v))*dx + p*div(v)*dx + q*div(u)*dx
L = inner(force, v) * dx \
+ pressure_left * inner(n, v) * ds(mark["left"]) \
+ pressure_right * inner(n, v) * ds(mark["right"])

noslip = Constant((0.0, 0.0))
bc_wall = DirichletBC(W.sub(0), noslip, subdomains, mark["wall"])

bc_left = DirichletBC(W.sub(1), pressure_left, subdomains, mark["left"])
bc_right = DirichletBC(W.sub(1), pressure_right, subdomains, mark["right"])

bcs = [bc_wall, bc_left, bc_right]

# Compute solution
w = Function(W)
solve(a == L, w, bcs)

# Split using deep copy
(u, p) = w.split(True)​


-Yang

Community: FEniCS Project

1 Answer


2
3 months ago by
There was a bug in FacetNormal implementation on quads in 2017.2.0 release. It is fixed in the development version and fix will appear in the upcoming release. Fix is here: https://bitbucket.org/fenics-project/ffc/commits/1473062346a43f591f3ec58ca92e572d2213a395.

(You can try a Docker container with the development version if it fixes your problem: https://quay.io/repository/fenicsproject/dev.)
Thank you so much for your help, Jan! That helps me a lot.
written 3 months ago by young zhou  
Hi Jan,

I followed your suggestion to try dev version out, but with quay.io/dolfinadjoint/dev-pyadjoint, I got the same unphysical result. While with quay.io/fenicsproject/dev, I got error "type object 'dolfin.cpp.mesh.CellType' has no attribute 'Type_quadrilateral'" appeared for 
mesh = UnitSquareMesh.create(30, 30, CellType.Type_quadrilateral)​
I aslo tried to generate mesh as listed below but it mentioned 'name 'UnitQuadMesh' is not defined'
mesh = UnitQuadMesh(mpi_comm_world(), 20, 20)​
 
written 3 months ago by young zhou  
In dev version:

UnitSquareMesh.create(32, 32, CellType.Type.quadrilateral)​

FEniCS team does not maintain pyadjoint images. You should contact dolfin-adjoint team.
written 3 months ago by Jan Blechta  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »