Hi Community,
I am trying to solve stokes equation with pure dirichlet bcs as follows:
Rectangle domain:
with zero velocity on top, bottom and around two circles, point(0.0, 0.0) has been pinned (for uniqueness problem) velocity on left and right equal to \boldsymbol{v}(x, y) = \big(\frac{6}{0.41^2}Uy(0.41 - y), 0\big) where U is velocity. Reynolds number for this problem = 1.0
# parameters
k = 2 #degree of polynomial
U = 15.0 #maximum velocity
mu = 1.0 #kinematic viscosity
circle1 = Circle (Point(0.2, 0.2), 0.05) #defining circle
circle2 = Circle (Point(1.2, 0.2), 0.05) #defining circle
Rec = Rectangle(Point(0.0, 0.0), Point(2.2, 0.41)) #defining rectangle
geometry = Rec - circle1 - circle2 #define geometry by subtracting rectangle and circle
mesh = generate_mesh(geometry, 80)
T = VectorElement("CG", mesh.ufl_cell(), 2)
Q = FiniteElement("CG", mesh.ufl_cell(), 1)
mixed = T * Q
W = FunctionSpace(mesh,mixed)
# defining subdomains
#defining pinpoint for removing degree of freedom for pressure
class PinPoint(SubDomain):
def inside(self, x, on_boundary):
return x[0] < DOLFIN_EPS and x[1] < DOLFIN_EPS
pinpoint = PinPoint()
# defining boundary condition around circle
inflow = 'near(x[0], 0)'
outflow = 'near(x[0], 2.2)'
topnbottom = 'near(x[1], 0) || near(x[1], 0.41)'
cylinder1 = 'on_boundary && x[0]>0.15 && x[0]<0.25 && x[1]>0.15 && x[1]<0.25'
cylinder2 = 'on_boundary && x[0]>1.15 && x[0]<1.25 && x[1]>0.15 && x[1]<0.25'
#No-slip boundary condition for velocity at top and bottom
noslip = Constant((0.0, 0.0))
bc0 = DirichletBC(W.sub(0), noslip, topnbottom)
#No-slip boundary condition for velocity around circle
bc1 = DirichletBC(W.sub(0), noslip, cylinder1)
bc2 = DirichletBC(W.sub(0), noslip, cylinder2)
#inflow bloundary condition for left and right boundary
infloww = Expression(('4.0*U*x[1]*(0.41 - x[1]) / pow(0.41, 2)', '0.0'), U = U, degree=4 + k)
bc3 = DirichletBC(W.sub(0), infloww, inflow)
zero = Constant((0.0))
bc4 = DirichletBC(W.sub(0), infloww, outflow)
bc5 = DirichletBC(W.sub(1), zero, pinpoint, "pointwise")
bcs = [bc0, bc1, bc2, bc3, bc4, bc5]
#Define variational problem
(v, p) = TrialFunctions(W)
(w, q) = TestFunctions(W)
f = Constant((0.0, 0.0))
a = mu*inner(grad(v), grad(w))*dx - div(w)*p*dx - q*div(v)*dx
L = inner(f, w)*dx
#Solver
V = Function(W)
solve(a == L, V, bcs)
#Get sub-functions
v, p = V.split()
#Save solution in VTK format
pfile_pvd = File("Steady_stokes/pressure.pvd")
pfile_pvd << p
I defined a test problem and got reasonable L2-error and convergence rate. but when I am gonna get pressure~~wrap by scalar I can’t see any symmetry in my result.
P.S. velocity looks fine in this case.
Thank you in advance.