How to apply pressure drop boundary conditions in Stokes flow?

General
1

Hi! I’m a new learner of FEniCS and trying to apply pressure drop boundary conditions to check if the Stokes flow is Darcy flow. I haven’t found any examples in the FEniCS demos. Please help me with where I can start. Any help is much appreciated. Here is my code.

import numpy as np
import time
from dolfin import ; from mshr import *
import dolfin.common.plotting as fenicsplot
from matplotlib import pyplot as plt
L = 4
H = 2
XMIN = 0.0; XMAX = L
YMIN = 0.0; YMAX = H
VIN = 4.0
resolution = 32
mesh = generate_mesh(Rectangle(Point(XMIN, YMIN), Point(L, H)) - Circle(Point(1.5, 0.25 * H), 0.2)
- Circle(Point(0.5, 0.5 * H), 0.2) - Circle(Point(2.0, 0.75 * H), 0.2)
- Circle(Point(3.0, 0.3 * H), 0.15) - Circle(Point(2.5, 0.6 * H), 0.1), resolution)
no_levels = 0
for i in range(0,no_levels):
cell_marker = MeshFunction(“bool”, mesh, mesh.topology().dim())
for cell in cells(mesh):
cell_marker[cell] = False
p = cell.midpoint()
if p.distance(Point(0.5, 0.5 * H)) < 0.5:
cell_marker[cell] = True
mesh = refine(mesh, cell_marker)
plt.figure()
plot(mesh, title=“Mesh”)
plt.show()
VE = VectorElement(“CG”, mesh.ufl_cell(), 2)
QE = FiniteElement(“CG”, mesh.ufl_cell(), 1)
WE = VE * QE
W = FunctionSpace(mesh, WE, constrained_domain=pbc)
V = FunctionSpace(mesh, VE, constrained_domain=pbc)
Q = FunctionSpace(mesh, QE, constrained_domain=pbc)
w = Function(W)
(u, p) = (as_vector((w[0],w[1])), w[2])
(v, q) = TestFunctions(W)
uin = Expression(("VIN4*(x[1](YMAX-x[1]))/(YMAXYMAX)", “0.”), YMAX=YMAX, VIN=VIN, element = V.ufl_element())
ib = Expression(“near(x[0],XMIN) ? 1. : 0.”, XMIN=XMIN, element = Q.ufl_element())
ob = Expression(“near(x[0],XMAX) ? 1. : 0.”, XMAX=XMAX, element = Q.ufl_element())
wb = Expression(“x[0] > XMIN + DOLFIN_EPS && x[0] < XMAX - DOLFIN_EPS ? 1. : 0.”, XMIN=XMIN, XMAX=XMAX, element = Q.ufl_element())
h = CellDiameter(mesh)
C = 1.0e10
gamma = C/h
f = Expression((“0.0”,“0.0”), element = V.ufl_element())
residual = ( - pdiv(v)dx + inner(grad(u), grad(v))dx + div(u)qdx +
gamma(ibinner(u - uin, v) + wbinner(u, v))*ds - inner(f, v)dx )
solve(residual == 0, w)
!rm results-NS/
plt.figure()
plot(u, title=“Velocity”)
plt.figure()
plot(p, title=“Pressure”)
plt.show()

2

Since you are new, make sure to read Read before posting: How do I get my question answered?. Furthermore, make sure to learn how to properly post your code: you must use “```”, and include all imports required by anyone else to be able to run your code. I removed your code since it didn’t comply with this rules, feel free to edit back the post with the corrected code.