Hi! I have to solve Stokes Equations on a cube domain (side L) with ellipsoidal obstacle in the middle of the cube (the cube and the obstacle have same center), with the following boundary conditions:
- null velocity on the surface of the obstacle
- null normal component of velocity on front and back face of the cube (z direction)
- mean of inflow pressure p1
- mean of outflow pressure null
- periodic condition for pressure and velocity on top and bottom face of the cube (y direction)
This is the code:
from dolfin import *
import scipy.io as io
import numpy as np
import math
#Set test variables
m = 1e-6;
ncells = 1;
eta = 0.001;
a = 7*m;
b = 10*m;
c = 10*m;
eps = 0.6649;
L = (4*math.pi*a*b*c/(3*(1-eps)))**(1.0/3.0);
#Inflow pressure
p1 = 1*ncells;
# Load mesh
mesh = Mesh("../Mesh/Ellipsoid.xml")
boundaries = MeshFunction("size_t", mesh, "../Mesh/Ellipsoid_facet_region.xml")
class PeriodicBoundary(SubDomain):
def inside(self, x, on_boundary):
return bool((x[1] < -L/2 + DOLFIN_EPS and
x[1] > (-L/2 - DOLFIN_EPS)) and on_boundary)
# Map right boundary (H) to left boundary (G)
def map(self, x, y):
y[0] = x[0]
y[1] = x[1] - L
y[2] = x[2]
# Define function spaces
V = VectorFunctionSpace (mesh, "CG", 2, constrained_domain=PeriodicBoundary())
Q = FunctionSpace (mesh, "CG", 1, constrained_domain=PeriodicBoundary())
W = MixedFunctionSpace([V,Q])
n = FacetNormal(mesh)
ds = Measure("ds", subdomain_data=boundaries)
# Create functions for boundary conditions
zero = Constant(0.0)
collector_bc = Constant((0,0,0))
# Boundary condition for velocity on extern boundaries
bc3 = DirichletBC(W.sub(0).sub(2), zero, boundaries, 3) #front
bc4 = DirichletBC(W.sub(0).sub(2), zero, boundaries, 4) #back
# Boundary condition for velocity on collector boundary
bc7 = DirichletBC(W.sub(0), collector_bc, boundaries, 7) #collector
# Collect boundary conditions
bcs = [bc3,bc4,bc7]
# Define variational problem
v, q = TestFunctions(W)
u, p = TrialFunctions(W)
a = eta*inner(grad(u), grad(v))*dx - div(v)*p*dx + q*div(u)*dx
L = - p1*dot(v,n)*ds(1)
# Compute solution
w = Function(W)
solve(a == L, w, bcs, solver_parameters={'linear_solver':'mumps'})
u, p = w.split(deepcopy=True)
# Save solution in VTK format
ufile_pvd = File("velocity.pvd")
ufile_pvd << u
pfile_pvd = File("pressure.pvd")
pfile_pvd << p
# Plot solution
plot(mesh)
plot(u)
plot(p)
interactive()
but the solution doesn't satisfy periodic boundary condition for velocity, as show in figure
And I don't find errors in my code.
Thank you in advance for your help.
Here the complete code with mesh.
