Relax a Neumann BC for Poisson problem


178
views
0
7 months ago by
Hi, I need to solve the Poisson problem with mixed BCs as shown below:

The red and black boundaries have natural (von Neumann) condition  $\nabla u\cdot\mathbf{n}=0$u·n=0 (black) and  $\nabla u\cdot n=1.0$u·n=1.0 (red). The red boundary also has the essential (Dirichlet) condition  $u=1.0$u=1.0.
I need to relax the essential boundary condition on the green parts, so that   $\nabla u\cdot n=unknown$u·n=unknown there.
The only restriction I can think of is that the total flow into the box balances the flow out of the box at the red part:$\int_{\partial\Omega_{green}}\nabla u\cdot n=-\int_{\partial\Omega_{red}}\nabla u\cdot n=-\int_{\partial\Omega_{red}}1.0_{ }$∂Ωgreenu·n=∂Ωredu·n=∂Ωred1.0

The physics here is potential flow (for ex https://en.wikipedia.org/wiki/Potential_flow) in a box with holes (green) and a moving piston on top (red).

I know how to add the condition for a known natural BC (as is done in demo_subdomains-poisson.py), but how do we let the gradient of  $u$u take any value on the green parts?
The code which does not let the gradient on the green parts vary is below, with black being the "Wall", red being the "Piston" and green being the "Holes".
Thank you!

from dolfin import *
import numpy as np

'''Solve poisson equation for potential flow in a 2d box with rising piston
'''

# Parameters for the boundary condition
y1 = 0.2
dhole = 0.1
dpiston = 0.2
vel_piston = 1.0
phi0piston = 1.0

class OnWall(SubDomain):
    def inside(self, x, on_boundary):
        # the boundary is everything but the side holes and top piston
        on_bottom = x[1] < DOLFIN_EPS + y1
        on_top = x[1] > (1.0 - DOLFIN_EPS - y1)
        on_side = x[1] < (1.0 - y1 - dhole) and x[1] > (y1 + dhole)
        return on_bottom or on_top or on_side

class OnHole(SubDomain):
    def inside(self, x, on_boundary):
        # the boundary is everything but the side holes and top piston
        on_bottom = x[1] < DOLFIN_EPS + y1
        on_top = x[1] > (1.0 - DOLFIN_EPS - y1)
        on_side = x[1] < (1.0 - y1 - dhole) and x[1] > (y1 + dhole)
        return (not (on_bottom or on_top or on_side))

class OnPiston(SubDomain):
    def inside(self, x, on_boundary):
        on_top = x[1] > (1.0 - DOLFIN_EPS - y1)
        return abs(x[0] - 0.5) < dpiston and on_top

# Create mesh and define function space
mesh = UnitSquareMesh(100, 100)
V = FunctionSpace(mesh, "Lagrange", 1)
u0 = Function(V)

# Initialize mesh function for boundary domains
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
wall = OnWall()
piston = OnPiston()
holes = OnHole()
wall.mark(boundaries, 0)
piston.mark(boundaries, 1)
holes.mark(boundaries, 3)
ds = Measure("ds")[boundaries]

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0.)
a = inner(nabla_grad(u), nabla_grad(v)) * dx
L = f * v * dx + vel_piston * v * ds(1)
# Define boundary conditions
# phi0top = Potential()
bc = DirichletBC(V, phi0piston, boundaries, 1)
# Compute solution
u = Function(V)
solve(a == L, u, bc)

# Hold plot
interactive()​

Community: FEniCS Project
Please login to add an answer/comment or follow this question.

Similar posts:
Search »