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

1 Answer


0
7 months ago by
Ovais  
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 \
        - dt * ke(Te, Tl) * inner(grad(Te), grad(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 \
        - dt * ke(Te, Tl) * inner(grad(Te), grad(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  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »