Using FEM to solve 1D coupled PDES [duplicate]

I have been attempting to numerically simulate the electron and lattice temperature of a two-layer metal subject to laser heating (1d two-temperature model). I have read related questions such as dealing with piece-wise coefficients and flux continuity; unfortunately, I'm still stumped with this code of mine as it involves coupled pdes that cover the entire region.

I'll just provide some context first to the problem: the main governing equations are given as follows:

$$C_e\frac{\partial T_e(z,t)}{\partial t}=\frac{\partial}{\partial z}\left(k_e\frac{\partial T_e(z,t)}{z} \right)- G(T_e(z,t)-T_l(z,t))+S(z,t)$$, $$C_l\frac{\partial T_l(z,t)}{\partial t}=\frac{\partial}{\partial z}\left(k_l\frac{\partial T_e(z,t)}{z} \right)+ G(T_l(z,t)-T_l(z,t))$$.

where $$T_e$$ and $$T_l$$ are the electron and lattice temperature of the entire multilayer, respectively; the functions I ultimately want. The parameters $$C,k,G,S$$ are the specific heat, thermal conductivity, coupling factor, laser source term, respectively (the subscripts are e="electrons" and l="lattice"). Because I want to investigate the heat transfer for 2 metals, I have all of these parameters to be position-dependent.

Here are the Neumann type boundary conditions: zero heat flux on the top ($$z=0$$) and bottom ($$z=b$$): $$\frac{\partial T_e}{\partial z}|_{z=0}=\frac{\partial T_l}{\partial z}|_{z=0}=0$$ and $$\frac{\partial T_e}{\partial z}|_{z=b}=\frac{\partial T_l}{\partial z}|_{z=b}=0.$$ The electron heat flux and temperature are also said to be continuous at the interface between the two metals. I assume FEM will be able to handle that part.

The initial conditions are $$T_e(z,t=0)=T_l(z,t=0)=300$$K.

Below is the simplified code that I have been working on. First, defining some laser parameters:

(*  required packages  *)

Needs["NDSolveFEM"]

(*  laser params *)

J = 30000.;(*laser fluence at 0 deg[J/m^2]*);
tp = 10*^-15;(*pulse duration[s]*)
\[Omega]0 = 10000.; (*laser beam waist*)
\[Lambda] = 512*^-9; (*wavelength*)
y = 0;
x = 0;
\[CurlyPhi] = 60*Pi/180; (*angle*)

d = 15.3*^-9;
Lambda]b = 100*^-9;

R = 0.93;

S[z_, t_] := 1/(tp (d + \[Lambda]b)) (0.94 (1 - Abs[R]^2) J Cos[\[CurlyPhi]]) Exp[-(z/d)] Exp[-2.77 (t/tp)^2] Exp[-2 (x^2/\[Omega]0^2 + y^2/(\[Omega]0^2/Cos[\[CurlyPhi]]))];


Here, I introduce the thermodynamic properties of each layer. Ideally, all of these would be functions of temperature/s, but I will set only Ce1 to be a function of $$T_e$$ and $$T_l$$:

(* layer 1 thermodynamic properties *)

Ce1 = (31.49683-33.33454/(1+Exp[(Te[t,z]/10000-0.87338)/0.31581]))\10^5;
ke1 = 2; (* random fillers for rest of params*)
Cl1 = 2;
kl1 = 2;
G1 = 2;

(* layer 2 thermodynamic properties *)

Ce2 = 95*((13.4*^11)*Te[t,z])/((7.9*^7)*(Te[t,z])^2+(13.4*^11)*Tl[t,z]);
ke2 = 3;
Cl2 = 3;
kl2 = 3;
G2 = 3;

(*layer params*)
a = 100*^-9; (* 100nm - first layer thickness *)
b = 100*^-9; (* 100nm - second layer thickness *)

tlimit = 15000*^-15;

T0 = 300; (* initial temp *)

(*property arrays*)
Ceset = {Ce1, Ce2};
keset = {ke1, ke2};
Clset = {Cl1, Cl2};
klset = {kl1, kl2};
Gset = {G1, G2};

(*Create mesh with region markers*)
g = {a, b} // N;(*thicknesses*)
gw = {0}~Join~Accumulate[g];
bmesh = ToBoundaryMesh["Coordinates" -> Partition[gw, 1],
"BoundaryElements" -> {PointElement[{{1}, {2}, {3}}]}]; nrEle = 10; \
pt = Partition[gw, 2, 1];
regmarkers =
Transpose[{Partition[(Mean /@ pt), 1], {1, 2},
Abs[Subtract @@@ pt]/nrEle}];
mesh = ToElementMesh[bmesh, "RegionMarker" -> regmarkers];

(*Create region dependent physical properties*)
Ce = Evaluate[
Piecewise[{{Ceset[[1]], ElementMarker == 1}, {Ceset[[2]],
ElementMarker == 2}}]];
ke = Evaluate[
Piecewise[{{keset[[1]], ElementMarker == 1}, {keset[[2]],
ElementMarker == 2}}]];
Cl = Evaluate[
Piecewise[{{Clset[[1]], ElementMarker == 1}, {Clset[[2]],
ElementMarker == 2}}]];
kl = Evaluate[
Piecewise[{{klset[[1]], ElementMarker == 1}, {klset[[2]],
ElementMarker == 2}}]];
G = Evaluate[
Piecewise[{{Gset[[1]], ElementMarker == 1}, {Gset[[2]],
ElementMarker == 2}}]];


And then here is the setup of the equations and conditions, which I am not sure if I have translated into code correctly; hence, no results (but also, no errors so I don't know what's wrong).

pde1 = D[ke*D[Te[t, z], z], z] - G*(Te[t, z] - Tl[t, z]) + Ss[t, z] -
Ce*D[Te[t, z], t];
pde2 = D[kl*D[Tl[t, z], z], z] + G*(Te[t, z] - Tl[t, z]) -
Cl*D[Tl[t, z], t];

ic = {Te[0, z] == Tl[0, z] == T0};

Subscript[\[CapitalGamma], N] =
NeumannValue[0, z == 0] + NeumannValue[0, z == Last@gw];

TTM2s = NDSolveValue[{pde1 == Subscript[\[CapitalGamma], N],
pde2 == Subscript[\[CapitalGamma], N], ic}, {Te[t, z],
Tl[t, z]}, {t, 0, tlimit}, {z} \[Element] mesh];

DensityPlot[Te[t, z], {t, 0, 15*^-12}, {z, First@gw, Last@gw},
ColorFunction -> "ThermometerColors", PlotRange -> All]


As mentioned, running the code didn't give any results nor errors. I'm not that well versed in Mathematica programming especially finite element methods. Can anyone identify what the faults are? Any help is appreciated!

Many thanks,

JR

• The output is inside TTM2s, not Te. Commented Nov 18, 2022 at 1:55
• Apologies @xzczd, a rookie mistake. I have a follow-up after trying to improve my code here. Commented Nov 19, 2022 at 15:10