### How to implement the two temperature model?

446

views

1

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!

\[

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

As a first step ,you can begin by understanding solution of coupled reaction diffusion equations. See https://fenicsproject.org/pub/tutorial/html/._ftut1010.html

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:

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.

r

`aise 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:

Happy new year!