How to implement the two temperature model?

446
views
1
7 months ago by
I would like to model heat diffusion at the gold / water interface after excitation of the metal surface by an ultrafast laser pulse (ca. 80 fs). I have already succeeded at implementing heat diffusion considering  only the rise of the lattice temperature, but a more appropriate model would be the "two temperature model" (TTM) or "parabolic two step" (PTS). This model considers the electronic and lattice systems as separate with electronic and lattice temperature $T_e$ and $T_l$, respectively. The electronic and lattice systems are coupled through an electron-phonon coupling constant $G$ and the electronic system is excited by a source term $S$ representing the incoming laser pulse. The model can be written as two coupled PDEs (from Smith (2001) Appl. Phys. Lett. 78:1240):
$C_e(T_e) \frac{\partial T_e}{\partial t} = \frac{\partial}{\partial x} \left( k_e(T_e, T_l) \frac{\partial T_e}{\partial x} \right) - G \left[ T_e - T_l \right] + S \\ C_l \frac{\partial T_l}{\partial t} = G \left[ T_e - T_l \right],$
where $C_e$, $C_l$ are electron and lattice heat capacities, respectively, and $k_e$ is the thermal conductivity with $k_e(T_e, T_l) = k_{eq}[T_e / T_l]$, $k_{eq} = 317$ W m$^{-3}$ K$^{-1}$ is the equilibrium thermal conductivity. The temperature-dependent electron heat capacity $C_e(T_e)$ may be approximated by $C_e(T_e) = C'_e T_e$ where $C'_e = 70$ J m$^{-3}$ K$^{-2}$.

Now, I am a novice, and I would like to ask how to go with the implementation of such kind of system where I have to consider the evolution of two coupled temperatures. Any code snippet, reference or reading material would be appreciated. Thanks in advance!
Community: FEniCS Project

0
7 months ago by
As a first step ,you can begin by understanding solution of coupled reaction diffusion equations. See https://fenicsproject.org/pub/tutorial/html/._ftut1010.html
Thanks for the reference! Using it, I tried to implement a MWE, but it fails with an arity mismatch:
raise ArityMismatch("Adding expressions with non-matching form arguments {0} vs {1}.".format(a, b))
ArityMismatch: Adding expressions with non-matching form arguments (Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2)), 1, None),) vs ().

Here's the MWE:
from fenics import *

# Coefficients
t_i = 0.               # initial time
t_f = 100e-12            # final time
time_zero = 1e-12
num_steps = 100      # number of time steps
dt = (t_f - t_i) / num_steps # time step size
initial_temperature = 300. # K

C_metal = .129      # heat capacity, gold / J g-1 K-1
rho_metal = 19.3    # density, gold / g cm-3

# Smith & Norris (2001) APL 78:1240
Ce_ = 70e-6            # J cm-3 K-2
keq = 317e-6           # W cm-3 K-1
G   = 2.1e10           # W cm-3 K-1
Cl  = C_metal * rho_metal

# Create mesh and define function space
nx = 16
ny = 16
width = 1e-3
height = 1e-3
mesh = RectangleMesh(Point(0., 0.), Point(width, height), nx, ny)

p = 1
V_ele = VectorElement("CG", mesh.ufl_cell(), p)
V = FunctionSpace(mesh, V_ele)

# Define boundary condition
u_D = Constant(initial_temperature)

tol = 1E-6 * height
# tol = 1E-12
class Boundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (near(x[0], 0, tol) or near(x[0], width, tol) or near(x[1], 0, tol) or near(x[1], height, tol))

boundary = Boundary()
boundary_markers = MeshFunction('size_t', mesh, 1)
boundary_markers.set_all(9999)
boundary.mark(boundary_markers, 0)

bcs = []
bcs.append(DirichletBC(V, Constant((300, 300)), boundary_markers, 0))

# We define two subdomains (metal == 0 and water == 1)
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
return x[1] <= height / 2. + tol

class Omega_1(SubDomain):
def inside(self, x, on_boundary):
return x[1] >= height / 2. - tol

# materials = CellFunction('size_t', mesh)
materials = MeshFunction('size_t', mesh, 2)    # NB: 'size_t' is the equivalent to 'uint' in newer Dolfin versions
subdomain_0 = Omega_0()
subdomain_1 = Omega_1()
subdomain_0.mark(materials, 0)
subdomain_1.mark(materials, 1)

# Defining the "trial" functions
(Te, Tl) = TrialFunctions(V) # these are Te_n+1 and Tl_n+1, respectively
(ve, vl) = TestFunctions(V)

# Define initial value
Te_n = interpolate(u_D, V.sub(0).collapse())    # these are Te_n and Tl_n, respectively
Tl_n = interpolate(u_D, V.sub(1).collapse())

# Define new measures
dx = Measure('dx', domain=mesh, subdomain_data=materials)
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)

Ce = lambda x: Ce_ * x
ke = lambda xe, xl: keq * xe / xl

F = (Ce(Te_n) * Te_n + dt * G * Tl) * ve * dx \
+ (Cl * Tl_n + dt * G * Tl) * vl * dx \
- Ce(Te) * Te * ve * dx \
- dt * G * Te * ve * dx \
- Cl * Tl * vl * dx \
- dt * G * Te * vl * dx

u = Function(V)
solve(F == 0, u, bcs)​

Happy new year!
written 7 months ago by François Lapointe
I found this question related to arity mismatch: https://www.allanswered.com/post/xngez/how-to-correctly-implement-integral-over-facets-in-dg/

I understand that the solution was to correctly separate the linear and bilinear forms. In my code, I lump everything into F and hope that the solver correctly chumps on it.

A potential solution would be to separate the linear from the bilinear forms beforehand:
F = (Ce(Te_n) * Te_n + dt * G * Tl) * ve * dx \
+ (Cl * Tl_n + dt * G * Tl) * vl * dx \
- Ce(Te) * Te * ve * dx \
- dt * G * Te * ve * dx \
- Cl * Tl * vl * dx \
- dt * G * Te * vl * dx

a = lhs(F)
A = assemble(a)

L = rhs(F)
b = assemble(L)

u = Function(V)
# solve(F == 0, u, bcs)
solve(A, u.vector(), b)​
However, this now fails with:
UFLException: Found Argument in denominator of <Division id=140653179287120> , this is an invalid expression.​
written 7 months ago by François Lapointe