Scalar wrapping of pressure for stokes problem

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.

Please take some time to carefully format your question such that code snippets are correctly presented.

Furthermore construct a minimal working example as described here.

Hi Dokken,

Thanks for your recommendation, I revised my question.

There are several things that are still missing, like the import statements and definition of k.

As it comes to symmetric solution, your problem is not symmetric? As you have y varying from 0 to 0.41, and the circles are centered at 0.2.

I have changed circle1 and circle2 definitions to:

cylinder1 = 'on_boundary && x[0]>0.155 && x[0]<0.255 && x[1]>0.155 && x[1]<0.255'
cylinder2 = 'on_boundary && x[0]>1.15 && x[0]<1.25 && x[1]>0.155 && x[1]<0.255'

but It seems that I am getting non-symmetric solution.

As I said in the previous post, as you havent made your example executable by missing definitions and import statements, I cannot investigate your problem. First thing you should do is to remove the circles, and see if you get reasonable pressure distribution in a normal channel flow. Then you should introduce one circle, and see if the results are reasonable.

You can also compare your solution to a modified version of the Stokes demo

A final way of verifying the implementation is the method of manufactured solutions, which is explained for a Stokes problem in this paper, page 17-18