Good morning everyone,
I am currently trying to model the cooling time for a square (so 2D) of a material and I am not very clear on how to implement the different subdomains.
To be specific: the system at hand should have Neuman BC (outward heat flux) on the perimeter of the square and Robin BC (outward convection) on the surface, within the perimeter (the latter excluded). Hopefully the picture will clarify (fluxes represented normal to the side/surface the exit):
This is my current implementation of the Neumann BCs:
L = 140E-6
#create the mesh, 2D
domain = Rectangle(Point(0.0, 0.0), Point(L, L))
mesh = generate_mesh(domain, 15)
V = FunctionSpace(mesh, 'P', 1)
#subdomain on which the flux is 0
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0, tol) or near(x[0], 1, tol)
#subdomain on which the flux is g
class Omega_1(SubDomain):
def inside(self,x,on_boundary):
return near(x[1], 0, tol) or near(x[1], 1, tol)
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1, 0)
subdomain_0 = Omega_0()
subdomain_0.mark(boundary_markers, 0)
subdomain_1 = Omega_1()
subdomain_1.mark(boundary_markers, 1)
#redefine boundary integration measure
ds = Measure('ds', domain = mesh, subdomain_data = boundary_markers)
then to actually apply the BC and compute the integrals over the boundary:
g = Constant(Q)
boundary_conditions = {0 : {'Neumann' : 0}, 1 : {'Neumann' : g}}
bcs=[] #empty BC list, no dirichlet
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(1/kappa*g*v*ds(i))
So, provided that what done until now is correct (on which I’d appreciate some feedback since I’m new with Fenics), how should I introduce a Robin BC that takes into account the convection from the square surface?
Thanks for your collaboration!