I’m sorry for this very specific problem but any help is appreciated.
I’m solving the time-harmonic heat equation
\begin{cases}
\kappa_\mathrm{e}\Delta T + i\omega\rho_\mathrm{e}c_\mathrm{e} T = 0 & \text{in } \mathcal{R}\setminus\overline{\mathcal{D}} \\
\kappa_\mathrm{i}\Delta T + i\omega\rho_\mathrm{i}c_\mathrm{i} T = 0 & \text{in }\mathcal{D} \\
T^+ = T^- & \text{on }\partial\mathcal{D}\\
-\kappa_\mathrm{e}\nabla T^+\cdot\mathbf{n} = -\kappa_\mathrm{i}\nabla T^-\cdot\mathbf{n} & \text{on }\partial\mathcal{D}\\
-\kappa_\mathrm{e}\nabla T\cdot\mathbf{n} = 0 & \text{on }\Gamma_\mathrm{sides}\\
-\kappa_\mathrm{e}\nabla T\cdot\mathbf{n} = h_\mathrm{eff}T- Q_\mathrm{lamp} & \text{on }\Gamma_\mathrm{front}\\
-\kappa_\mathrm{e}\nabla T\cdot\mathbf{n} = h_\mathrm{eff}T & \text{on }\Gamma_\mathrm{back}\\
\end{cases}
in a plate
\mathcal{R}=\left\{\left(x,y\right)\in\mathbb{R}^2\,:\,0<x<L_x,\,0<y<L_y\right\}
(L_x=0.01, L_y=1) which is irradiated on one of its sides:
\Gamma_{\mathrm{front}}=\left\{ (0,y)\,\vert\; 0 < y < L_y \right\}
The function \mathcal{Q}_\mathrm{lamp} is:
\mathcal{Q}_\mathrm{lamp} =\alpha\frac{ I}{2\pi}\frac{\mathbf{x}\cdot\mathbf{i}}{\left(\mathbf{x}-\mathbf{x}_l\right)\cdot\left(\mathbf{x}-\mathbf{x}_l\right)}
being \mathbf{x}_l the location from where the plate is being irradiated:
For solving that problem I made the following code (which takes the mesh which was generate by gmsh):
def Thermogram(mesh_and_tags, x_l, y_l, w, I=1000 ):
'''Thermogram(mesh_and_tags, x_l, y_l, w ) -> T
where
mesh_and_tags = (mesh, cell_tags, facet_tags)
being the tags:
1 : plate
2 : front
3 : back
4 : sides
5 : defect
'''
mesh, cell_tags, facet_tags = mesh_and_tags
ds = Measure("ds", domain=mesh, subdomain_data=facet_tags)
tags = { "plate" : 1,
"front": 2,
"back" : 3,
"sides" : 4,
"defect" : 5}
# Physical parameters:
alpha = 0.3
h = 15
To = 298
h_eff = h + alpha*sigma*4*To**3
k_e = 200 #
k_i = 0.025 #
c_e = 900
c_i = 1000
rho_e = 2700
rho_i = 1
C0 = FunctionSpace(mesh, ("DG", 0))
kappa = Function(C0)
iwrhoc = Function(C0)
kappa.x.array[ cell_tags.find(tags["plate"])] = np.full_like(cell_tags.find(tags["plate"]), k_e, dtype=ScalarType)
kappa.x.array[ cell_tags.find(tags["defect"])] = np.full_like(cell_tags.find(tags["defect"]), k_i, dtype=ScalarType)
iwrhoc.x.array[ cell_tags.find(tags["plate"])] = np.full_like(cell_tags.find(tags["plate"]), 1j*w*rho_e*c_e, dtype=ScalarType)
iwrhoc.x.array[ cell_tags.find(tags["defect"])] = np.full_like(cell_tags.find(tags["defect"]), 1j*w*rho_i*c_i, dtype=ScalarType)
H = FunctionSpace(mesh,('Lagrange',1))
phi = TrialFunction(H)
psi = TestFunction(H)
T = Function(H)
T.name = 'T'
Q = Function(H)
Q.name = 'Q'
Q.interpolate( lambda x : alpha*I/(2*np.pi)*( x[0] - x_l )/ ( (x[0]-x_l)**2 + (x[1] - y_l)**2 ))
# variational formulation
a = kappa*inner(grad(phi), grad(psi))*dx + iwrhoc*inner(phi, psi)*dx + h_eff*inner(phi,psi)*ds(tags["front"]) + h_eff*inner(phi,psi)*ds(tags["back"])
l = h_eff*inner(Q,psi)*ds(tags["front"])
problem = LinearProblem(a, l)
T = problem.solve()
return T, H
My problem is that, when testing (for example with a lamp at (x,y) = (-0.2,0.2)) I get funny behaviour .
For example, at x=0 I get:
The spike is expected due to the change in the thermal properties in \mathcal{D}, however at the left, it looks as if the second derivative wasn’t monotonic. Indeed, it doesn’t seem to be fulfilling the homogeneous Neumann boundary conditions con \Gamma_{\mathrm{sides}}
If you need the code for the mesh generation I can add it too.
EDIT: I add the code for generating the mesh:
from numpy import pi
import pygmsh
def make_real_domain( xy_defect, r, Lx = 0.01, Ly = 1, Nx = 5, Nd = 30):
'''
make_real_domain( xy_defect, r, Lx = 0.01, Ly = 1, Nx = 5, Nd = 30)
xy_defect are given as the fraction with respect to Lx and Ly
r is given as a fraction of Lx
returns:
mesh, cell_tags, facet_tags
being the tags:
1 : plate
2 : front
3 : back
4 : sides
5 : defect
'''
res_plate = Lx / Nx
# defect parameters
x_d , y_d = xy_defect
x_d = x_d*Lx
y_d = y_d*Ly
r = r*Lx
res_circ = 2*pi*r / Nd
# Initialize empty geometry using the build in kernel in GMSH
geometry = pygmsh.geo.Geometry()
# Fetch model we would like to add data to
model = geometry.__enter__()
# Add defect
defect = model.add_circle( (x_d, y_d), r, mesh_size=res_circ)
# Add corners of the plate
points = [model.add_point(( 0, 0, 0), mesh_size=res_plate),
model.add_point((Lx, 0, 0), mesh_size=res_plate),
model.add_point((Lx, Ly, 0), mesh_size=res_plate),
model.add_point(( 0, Ly, 0), mesh_size=res_plate)]
# Add lines between all points creating the rectangle
plate_boundaries = [model.add_line(points[i], points[i+1])
for i in range(-1, len(points)-1)]
# Create a line loop and plane surface for meshing
plate_boundary = model.add_curve_loop(plate_boundaries)
plate_surface = model.add_plane_surface(plate_boundary, holes=[defect.curve_loop])
def_surface = model.add_plane_surface(defect.curve_loop)
# Call gmsh kernel before add physical entities
model.synchronize()
# not sure if needed
model.add_physical([plate_surface], "Plate")
model.add_physical([plate_boundaries[0]], "Front")
model.add_physical([plate_boundaries[2]], "Back")
model.add_physical([plate_boundaries[1], plate_boundaries[3]], "Sides")
model.add_physical([def_surface], "Defect")
dim = 2
geometry.generate_mesh(dim=dim)
filename = 'temp.msh'
import gmsh
gmsh.write(filename)
gmsh.clear()
geometry.__exit__()
from dolfinx.io import gmshio
from mpi4py import MPI
gmsh.initialize()
gmsh.model.add("Mesh from file")
gmsh.merge(filename)
mesh, cell_tags, facet_tags = gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=dim)
gmsh.finalize()
return ( mesh, cell_tags, facet_tags)