Choosing a right approach to a problem (heat equation + Ohm's law in 2 adjacent regions)

I would like to simulate a physical experiment known as the 3 omega method which consists of using an AC current through a wire deposited onto the surface of a thin film. Then, the voltage of the wire is read, and should give information about the thermal conductivity of the thin film, and this is what I would like to see with FEniCS.

I have created a mesh with Gmsh that consists of 2 parallelepipedes (2 boxes using the OpenCascade environment), one on top of the other. The thin film is a thin flat rectangular box, while the wire is a long thin square cross section box. I have used boolean fragments to ensure that the surface they have in common is meshed only once. All seems to work fine regarding this aspect.

Then here comes the doubt. My approach would be to use a mixed element space, at least for the region of the wire, where I have to solve a heat equation that involves a Joule term. I will have to assume a temperature dependent (linear) electric resistance, so it is question to solve the heat equation containing a Joule term and also to solve for the voltage, thus I must also solve Ohm’s law in the material, which is another weak formulation problem, hence the mixed element approach (is this approach justified?).
For the thin film domain, I plan to only solve for the heat equation that contains only the transient term (involving the specific heat of the thin film) and the Fourier conduction term (involving its thermal conductivity), as almost no electric current should flow through the film, hence Joule heat is negligible and there is essentially no voltage through the film. The film is just in thermal contact with the wire. So, for the thin film region, I would consider its bottom surface to be at a fixed temperature (Dirichlet boundary condition).

Hence the weak form for the heat equation of the whole problem would contain terms that involve, say dx1 for the wire region, and dx2 for the thin film region. The only boundary condition would be Dirichlet at the bottom surface of the thin film (i.e. region 2). Ohm’s law would be solved with Neumann or Dirichlet boundary conditions to make sure a specified current is passing through the wire, and this would be solved only in region 1, so it would only involve dx1 but no dx2.

Does this seem like a correct approach?

Edit:

I have faced some difficulties and while I am almost there, I am still facing a problem: Newton method does not converge and I think it isn’t because of the initial guess… but because the code is not doing what I think it does, so a part is likely wrong.

I have learned the hard way that I cannot follow the approach I describe above. It seems that if I solve a heat equation in a domain, then I must also solve a heat equation in an adjacent domain. It is disappointing to me because I would have liked to solve the PDE only in a single domain (Ohm’s law actually), since the electric current is so small in the thin film that it shouldn’t have a significant impact on the solution anyway. That is because the metal wire conducts electricity much better than the thin film, so most of the current passes through the wire, and very little current extends into the thin film region. So I am not even interested to solve for Ohm’s law (to retrieve the voltage for example, or the current) inside the thin film, but if I don’t, I get an error.

So here is my full code, currently it fails because Newton method fails. I am running out of ideas to debug the code, i.e. find the term(s) that are problematic/wrong. What approach would you recommend to debug my code?

import numpy as np
import meshio
from scipy import interpolate as sci_interpolate

msh = meshio.read("three_omega_thin_film.msh")

cells = np.vstack(np.array([cells.data for cells in msh.cells if cells.type == "tetra"]))
tetra_data = msh.cell_data_dict["gmsh:physical"]["tetra"]
tetra_mesh = meshio.Mesh(points=msh.points, cells=[("tetra", cells)], cell_data={"name_to_read": [tetra_data]})
facet_cells = np.vstack([cells.data for cells in msh.cells if cells.type == "triangle"])
facet_data = msh.cell_data_dict["gmsh:physical"]["triangle"]
facet_mesh = meshio.Mesh(points=msh.points,
                         cells=[("triangle", facet_cells)],
                         cell_data={"name_to_read": [facet_data]})
# Write mesh
meshio.xdmf.write("mesh_test0.xdmf", tetra_mesh)
meshio.xdmf.write("facet_mesh_test0.xdmf", facet_mesh)

# Use dolfin below.
from dolfin import *
mesh = Mesh()

# Read mesh.
with XDMFFile("mesh_test0.xdmf") as infile:
    infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("facet_mesh_test0.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
mvc2 = MeshValueCollection("size_t", mesh, 3)
with XDMFFile("mesh_test0.xdmf") as infile:
    infile.read(mvc2, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc2)
dx = Measure("dx", domain=mesh,subdomain_data=cf)
# External surface
ds = Measure("ds", domain=mesh, subdomain_data=mf)
# Internal surface
dS = Measure("dS", domain=mesh, subdomain_data=mf)
subdomains = cf
boundary_parts = mf

# Volumes
# 1 thin film
# 2 heater wire

# Surfaces
# 3 left end heater wire.
# 4 right end heater wire.
# 5 Interface between wire and thin film.
# 6 Thin film bottom


dx_thin_film = dx(subdomain_id=1)
dx_wire = dx(subdomain_id=2)
ds_bottom_thin_film = ds(subdomain_id=6)
dS_interface=dS(subdomain_id=5)

print('Expected thin film volume:', 1e-2 * 5e-3 * 1e-4)
print('Computed volume: ',assemble(Constant(1)*dx_thin_film))
#print('computed volume again: ',assemble(Constant(1)*dx(domain=mesh, subdomain_data=cf, subdomain_id=1)))

print('Expected wire volume:', 5e-3 * 1e-4 * 1e-4)
print('Computed volume: ',assemble(Constant(1)*dx_wire))

print('Expected interface area:', 5e-3* 1e-4) 
print('Computed surface area: ',assemble(Constant(1)*dS_interface))

print('Expected bottom thin film surface: ', 1e-2 * 5e-3)
print('Computed bottom thin film surface: ', assemble(Constant(1)*ds_bottom_thin_film))

P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
# Mixed function space
ME = FunctionSpace(mesh, P1*P1)

# Physical parameters
T_ref = 300.0
kappa_wire = Constant(75.0)
kappa_thin_film = Constant(7.8)
rho_interpolated = sci_interpolate.InterpolatedUnivariateSpline([10,400], [2.25e-9, 1.14e-7], k=1) # linear interpolation
current = Constant(90e-1)
J_0 = current / assemble(1 * ds(3))

# Steady-state solution first, to be used as initial guess of the transient solution.
dTsteadyVsteady = TrialFunction(ME)
usteady, vsteady = TestFunctions(ME)            
TsteadyVsteady = Function(ME)
Tsteady, Vsteady = split(TsteadyVsteady)

V = FunctionSpace(mesh, 'CG', 1)
rho_steady = project(float(rho_interpolated(T_ref)),V)
rho_thin_film = project(5.0e-3,V)

TsteadyVsteady_init=Expression((T_ref,'0.225'),degree=1)
TsteadyVsteady.interpolate(TsteadyVsteady_init)
J_injected = J_0
J_local = -(1.0/rho_steady)*grad(Vsteady)
J_thin_film = -(1.0/rho_thin_film)*grad(Vsteady)

Fourier_terms = -dot(kappa_wire * grad(Tsteady), grad(usteady)) * dx_wire - dot(kappa_thin_film * grad(Tsteady), grad(usteady)) * dx_thin_film
Joule_term = dot(rho_steady*J_local, J_local) * usteady*dx_wire + dot(rho_thin_film*J_thin_film, J_thin_film) * usteady*dx_thin_film# Add Joule heat in thin film.

FsteadyT = Fourier_terms + Joule_term
FsteadyV = dot(rho_steady * grad(Vsteady), grad(vsteady))*dx_wire + J_injected * vsteady * ds(4) + dot(rho_thin_film * grad(Vsteady), grad(vsteady))*dx_thin_film

weak_form = FsteadyT + FsteadyV
Jacsteady = derivative(weak_form,TsteadyVsteady,dTsteadyVsteady)

V_left = Constant(0.0)
bc0 = DirichletBC(ME.sub(0), T_ref, boundary_parts, 6)
bc1 = DirichletBC(ME.sub(1), V_left, boundary_parts, 3)
bcs = [bc0, bc1]

problem = NonlinearVariationalProblem(weak_form, TsteadyVsteady, bcs, J = Jacsteady)
solver = NonlinearVariationalSolver(problem)
solver.solve()

The file.msh is created with the following code and gmsh:

//+
SetFactory("OpenCASCADE");
Box(1) = {0, 0, 0, 1e-2, 5e-3, 1e-4};
//+
Box(2) = {2.5e-3, 2.5e-3, 1e-4, 5e-3, 1e-4, 1e-4};

//+
BooleanFragments{ Volume{1}; Delete; }{ Volume{2}; Delete; }
//+
Physical Volume("thin_film") = {1};
//+
Physical Volume("heater_wire") = {2};
//+
Physical Surface("heater_end_left") = {7};
//+
Physical Surface("heater_end_right") = {8};
//+
Physical Surface("wire_and_thin_film_interface") = {11};

//+
MeshSize {19, 23, 21, 17, 18, 20, 24, 22} = 5e-4;
//+
MeshSize {11, 9, 10, 12, 15, 13, 16, 14} = 1e-4;
//+
Physical Surface("thin_film_bottom") = {17};

Edit 2: If I comment out the Joule term, the code “works”, although obviously the wire doesn’t heat, so the temperature is homogeneous. And this seems to point at a clue, because with Paraview I can see that the voltage seems erroneous too: it is almost zero everywhere but at one end of the wire, where it is worth -7e11, a gigantic value.

Edit: As pointed out in one of the replies, I spotted a mistake in the weak form for the voltage, there should be a 1/rho instead of rho. Now the code converges towards realistic values! I still need to work on it and make rho a function of temperature rather than a constant, but otherwise all looks fine to me.

Have you considered the Mixed FEM approach from : https://arxiv.org/pdf/1911.01166.pdf and the examples within ?

No, I wasn’t aware of this paper, thanks for the reference. So, you think the problem might be the way I employ the mixed elements?

Yes I do, with mesh views you can solve equations on parts of the domain and couple them.

I see, thanks for the info. I will keep this in mind if I really need a fast execution code, in which case I might go through the hassle to get the code working to solve the PDE in a single region.
But as you can see in my edit to the OP, I have decided to just solve the two PDEs in both regions. My problem is that Newton method is failing for an unknown reason to me.

As an attempt to debug, I have also simplified my code and focused solely on the heat equation with the Joule terms (so no mixed elements), and it works as expected. So I know the mesh and its regions are fine.

I think I have spotted something wrong. The way I defined rho_steady, as rho_steady = project(float(rho_interpolated(T_ref)),V), i.e. I use a projection over “V”, a function space not related in any way to the ME. I have tried to replace “V” by “ME.sub(1)” for instance, but this returns an error.
Do you have an idea on how I could fix this problem? If I pick a constant instead (rho = Constant(1e-4)) then the code is able to run, but the values returned are totally off, I believe. It’s as if only the right end side of the wire is heated up by Joule heat, instead of the whole wire.

Edit: I have spotted a mistake regarding the weak form for the voltage equation. It should be 1/rho rather than rho. Now the code converges towards realistic values! I just need to make sure that rho is a function of T rather than a constant, but otherwise all looks OK to me!