Scalar wrapping of pressure for stokes problem

1

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.

pressure
pressure1362×860 52.2 KB

P.S. velocity looks fine in this case.

Thank you in advance.

2

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.

3

Hi Dokken,

Thanks for your recommendation, I revised my question.

4

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.

5

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.

pressure_n
pressure_n966×611 88.9 KB

6

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