Robin condition on complex shape

1

Hello, I am trying to study heat diffusion in a radiator composed of fins.
I began by creating a basic mesh that I can duplicate to create multiple fins, then added a Neumann condition on the bottom and a Robin condition anywhere else.
This is my code

def generateFin(z,lx,ly,lp,e):
    support = Box(Point(lx/2-lp/2,ly/2-lp/2,z),Point(lx/2+lp/2,ly/2+lp/2,z+e))
    fin = Box(Point(0.0,0.0,z+e),Point(lx,ly,z+2*e))
    return fin+support

fin = generateFin(0,lx,ly,lp,e)
fin2 = generateFin(2*e,lx,ly,lp,e)
mesh = generate_mesh(fin, 20)

plt.figure(1)
plot(mesh, title="Mesh")


Hh2 = FunctionSpace(mesh, 'P', 3)


#Boundary conditions

boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)

class BoundaryNeumann(SubDomain): 
    def inside(self, x, on_boundary):
        return on_boundary and near(x[2], 0, tol)
S = BoundaryNeumann()
S.mark(boundary_markers, 0)


class BoundaryRobin(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and x[2] > tol
    
S = BoundaryRobin()
S.mark(boundary_markers, 1)


boundary_conditions = {0: {'Neumann': j_a},
                       1: {'Robin': (h_a, thetaA)}}


ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)

u = TrialFunction(Hh2)
v = TestFunction(Hh2)

bcs = []
for i in boundary_conditions:
    if 'Dirichlet' in boundary_conditions[i]:
        bc = DirichletBC(Hh2, boundary_conditions[i]['Dirichlet'],
                         boundary_markers, i)
        bcs.append(bc)

integrals_N = []
for i in boundary_conditions:
    if 'Neumann' in boundary_conditions[i]:
        if boundary_conditions[i]['Neumann'] != 0:
            g = boundary_conditions[i]['Neumann']
            integrals_N.append(g*v*ds(i))

integrals_R_a = []
integrals_R_L = []
for i in boundary_conditions:
    if 'Robin' in boundary_conditions[i]:
        r, s = boundary_conditions[i]['Robin']
        integrals_R_a.append(r*u*v*ds(i))
        integrals_R_L.append(r*s*v*ds(i))


A = Expression(((str(k), '0', '0'), ('0', str(k), '0'), ('0','0',str(k))), degree=1)
a = dot(dot(A,grad(u)), grad(v))*dx + sum(integrals_R_a)
l = sum(integrals_N) + sum(integrals_R_L)

u = Function(Hh2)
solve(a == l, u)

It basically works (results below, I can only upload one image)

Then I try to add a second fin, by editing
mesh = generate_mesh(fin, 20)
to
mesh = generate_mesh(fin+fin2, 20)

But now, I get these very strange results, my mesh looks good but I suspect my boundary conditions not to be correctly defined, even though it worked fined with only one fin. What did I do wrong ? Is there another way to achieve what I want to do ?

image
image883×484 120 KB

2

As you have not supplied a complete code (There are several import statements and definitions of variables missing), It is tricky to give you accurate advice.
Please supply a Minimal Working Example, that you make sure can be run by others.

However, I can give you some pointers:

  • To check that you have marked your boundaries correctly, you can visualize the boundary markers, (for instance with Paraview, by saving them to pvd. File("mf.pvd").write(boundary_markers)

  • Additionally, your marking of Neumann boundaries has no effect, as you have initialized the mesh function with a 0 value everywhere. I would suggest to change the markings to use 1 and 2 for Robin and Neumann conditions respectively.

3

Thanks for your answer, here is my complete code:

from dolfin import *
from fenics import *
from mshr import *
import matplotlib.pyplot as plt

#Dimensions
lx=0.15
ly=0.08
lz=0.15
e = 0.01
lp = 0.038
#Physics
lbda=220
c=900
rho=2700
flux = 120/10
j=flux/(lp**2) 
h=150
Tobs = 1
theta0 = 1

k = lbda*Tobs/(rho*c)
j_a = j*Tobs*theta0/(rho*c)
h_a = h*Tobs/(rho*c)
thetaA=303.15
tol = 1E-14


###Mesh creation and vizualisation###

def generateFin(z,lx,ly,lp,e):
    support = Box(Point(lx/2-lp/2,ly/2-lp/2,z),Point(lx/2+lp/2,ly/2+lp/2,z+e))
    fin = Box(Point(0.0,0.0,z+e),Point(lx,ly,z+2*e))
    return fin+support

fin = generateFin(0,lx,ly,lp,e)
fin2 = generateFin(2*e,lx,ly,lp,e)
mesh = generate_mesh(fin+fin2, 20)

plt.figure(1)
plot(mesh, title="Mesh")


###Function space###
Hh2 = FunctionSpace(mesh, 'P', 3)


###Boundary conditions###
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)

#Bottom surface

class BoundaryNeumann(SubDomain): 
    def inside(self, x, on_boundary):
        return on_boundary and near(x[2], 0, tol)
S = BoundaryNeumann()
S.mark(boundary_markers, 1)

#Anywhere else

class BoundaryRobin(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and x[2]>tol
    
S = BoundaryRobin()
S.mark(boundary_markers, 2)


boundary_conditions = {1: {'Neumann': j_a},
                       2: {'Robin': (h_a, thetaA)}}


###variational formulation###

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)

u = TrialFunction(Hh2)
v = TestFunction(Hh2)

integrals_N = []
for i in boundary_conditions:
    if 'Neumann' in boundary_conditions[i]:
        if boundary_conditions[i]['Neumann'] != 0:
            g = boundary_conditions[i]['Neumann']
            integrals_N.append(g*v*ds(i))

integrals_R_a = []
integrals_R_L = []
for i in boundary_conditions:
    if 'Robin' in boundary_conditions[i]:
        r, s = boundary_conditions[i]['Robin']
        integrals_R_a.append(r*u*v*ds(i))
        integrals_R_L.append(r*s*v*ds(i))


A = Expression(((str(k), '0', '0'), ('0', str(k), '0'), ('0','0',str(k))), degree=1)
a = dot(dot(A,grad(u)), grad(v))*dx + sum(integrals_R_a)
l = sum(integrals_N) + sum(integrals_R_L)

u = Function(Hh2)
solve(a == l, u)

###Matplotlib for quick vizualisation###
'''
plt.figure(2)
plu=plot(u-273.15)
plt.title('Output')
plt.xlabel('x')
plt.ylabel('y')
plt.colorbar(plu)
'''
###File output###
'''
File("mf.pvd").write(boundary_markers)
file = File ("output2.pvd")
file << u
'''

I tried to view my boudaries markers with paraview, it seems correct. I also tried to change the markings as you suggested without success.

EDIT: I think I found what was happening, I used a too weak resolution in my generate_mesh function, which caused border conditions not to be properly defined, it now works with an higher resolution.