2D heat equation with 2 Neumann & 2 Diri

1

I try to run the code below (2D steady state heat equation) with 2 Neumann conditions on left and right respectively and 2 dirichlet on top and bottom respectively. I followed a previous post explaining how to implement multiple Neumann Boundaries (dirichlet boundary conditions are relatively easy to do).

However the code below does not work properly: I succeed to apply the Neumann condition on Left but not on right.

Any idea to find the error ?
Also if I want to apply a different flux on each of the two Neumann Boudaries, is my code correctly made ?

Thanks

#Plate with Boundary Conditions (BC): psi (Flux) on left & right sides, theta_0 (temperature) on top & Bottom sides
#C1 Load modules
from dolfin import *
from future import print_function
from fenics import *
import matplotlib.pyplot as plt

kappa = 1 #non dimensional diffusivity
psi = 10 #non dimensional flux
theta_0 = 10. # Celcius, temperature
tol = 1E-14

#C3 Geometry of the problem
X=1
Y=1

#C4 Mesh
nx2=20
ny2=20

mesh2 = RectangleMesh(Point(0.0, 0.0), Point(X, Y), nx2, ny2, diagonal=“left”)

plt.figure(1)
plot(mesh2, title=“Mesh of Plate”)

#C5 Space definition to search the solution
Hh2 = FunctionSpace(mesh2, ‘P’, 1)

#C6 Boundary Conditions

#C6.1 Top side BC: Dirichlet BC → Temperature imposed
def Top(x, on_boundary):
return on_boundary and (abs(x[1] - 1) < tol)
reD2 = DirichletBC(Hh2, Constant(theta_0), Top)

#C6.1bis Bottom side BC: Dirichlet BC → Temperature imposed
def Bot(x, on_boundary):
return on_boundary and (abs(x[1] - 0) < tol)
reD4 = DirichletBC(Hh2, Constant(theta_0), Bot)

bc = [ reD2, reD4 ]

######################################################
# Neumann BC on left, right and bottom sides → Flux imposed
# Create classes for defining parts of the boundaries
class Left(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (abs(x[0] - 0) < tol)

class Right(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (abs(x[0] - 1) < tol)

    # Initialize sub-domain instances

leN2 = Left()
RiN2 = Right()

    # Initialize mesh function for boundary domains

boundaries2 = MeshFunction(“size_t”, mesh2, mesh2.topology().dim()-1, 0 )
boundaries2.set_all(0)
leN2.mark(boundaries2, 1)
RiN2.mark(boundaries2, 2)

    # Define new measures associated with the exterior boundaries

ds2 = Measure(“ds”, domain=mesh2, subdomain_data=boundaries2)
#Fin C6.2 Left side: Neumann BC → Flux imposed
#############################################################

C7 Variationnal Formulation

u2 = TrialFunction(Hh2)
v2 = TestFunction(Hh2)
l2 = psiv2ds2(1)
a2 = kappa*inner(grad(u2), grad(v2))*dx

C8 Find the solution (SOLVER)

u2 = Function(Hh2)
solve(a2 == l2, u2, bc)

C8 Plot the solution

plt.figure(2)
plu=plot(u2)
plt.title(‘Numerical Solution u2’)
plt.xlabel(‘x’)
plt.ylabel(‘y’)
plt.colorbar(plu)

2

Please format your code appropriately using 3x` encapsulation, so your code is formatted properly.

There are also plenty of examples on this forum on how to debug your code when it comes to boundary conditions. First of all your MeshFunction (boundaries2) to file and visualize it with Paraview to make sure that the correct boundaries are marked.