How to implemenent Interface Resistance between metal non-metal layer

[Sandip433]

Hi, I want to simulate effect of heat transfer across Metal-Nonmetal Interface on femtosecond Laser. Here are the equations I want to solve
Interface Resistance.pdf
Here is the code I wrote

from fipy import *

Lx=2.0e-7 
dx=1.0e-10
nx=int(Lx/dx) 
T=100e-12
dt=4e-15
I_0 = 1100 
R = 0.982
delta = 12.23e-9
t_p = 100e-15
m=Grid1D(nx=nx, dx=dx)
x = m.cellCenters[0]
time = Variable(0.)
gold = (x <= Lx / 4.0)
chromium = (x >= Lx / 4.0 ) & ( x<= Lx/2.0)
metal = (x<= Lx/2.0)
nonmetal = (x>= Lx/2.0)
S = numerix.sqrt(4*numerix.log(2)/numerix.pi)*((1-R)*I_0/(t_p*delta))*numerix.exp(-(x/delta)-4*numerix.log(2)*((time-2*t_p)/t_p)**2)
T_e=CellVariable(name="Electron Temperature", mesh=m, value=300.0, hasOld=True )
T_p=CellVariable(name="Phonon Temperature", mesh=m, value=300.0, hasOld=True )
T_s=CellVariable(name="Non-metal Temperature", mesh=m, value=300.0, hasOld=True )
C_e = (68 * gold + 194 * chromium + 0 * nonmetal) * T_e
C_p = (2.5e6 * gold + 3.3e6 * chromium + 0 * nonmetal)
C_s = (0 * gold + 0 * chromium + 2e6 * nonmetal)
G_g = 0.21e17*((1.18e7/1.25e11)*(T_e+T_p)+1)
G_c = 4.2e17*((7.9e7/13.4e11)*(T_e+T_p)+1)
K_eg = 318*(1.25e11*T_e/(1.18e7*T_e**2+1.25e11*T_p))
K_ec = 95*(13.4e11*T_e/(7.9e7*T_e**2+13.4e11*T_p))
K_e = K_eg * gold + K_ec * chromium + 0 * nonmetal
K_s = (0 * gold + 0 * chromium + 10 * nonmetal)
G = G_g * gold + G_c * chromium + 0 * nonmetal
M=1e10 # taking R_es=R_ps=1e10
eq_e=TransientTerm ( coeff= C_e , var=T_e) == DiffusionTerm (coeff=K_e, var=T_e) - ImplicitSourceTerm (coeff=G, var=T_e) + ImplicitSourceTerm (coeff=G, var=T_p) + S * gold
eq_p=TransientTerm ( coeff= C_p , var=T_p) == DiffusionTerm (coeff=(0.01 * K_e ) , var=T_p) + ImplicitSourceTerm (coeff=G, var=T_e) - ImplicitSourceTerm (coeff=G, var=T_p)
eq_s=TransientTerm ( coeff= C_s , var=T_s) == DiffusionTerm (coeff=K_s, var=T_s) + ImplicitSourceTerm (coeff=  M * nonmetal, var=T_e) - ImplicitSourceTerm (coeff= M * nonmetal, var=T_s) + ImplicitSourceTerm (coeff= M * nonmetal, var=T_p) - ImplicitSourceTerm (coeff=M * nonmetal, var=T_s) 
eq=eq_e & eq_p & eq_s
solver = LinearPCGSolver(precon=None, tolerance=1e-10, iterations=2000)
while time()<T:  
    T_e.updateOld()
    T_p.updateOld()
    T_s.updateOld()
    for i in range(10):
        res=eq.sweep(dt=dt, solver=solver)
    time.setValue(time() + dt)

The problem is when I solve those equations for 100 ps the temperature of non-metal substrate does not change. What should I do ?

[guyer]

Why would T_s ever change? The equation you have written is (including source terms that are not in your PDF)

$$\frac{\partial C_s T_s}{\partial t} = \nabla\cdot\left(k_s\nabla T_s\right) + M \cdot {\tt nonmetal} \left(T_e - T_s + T_p - T_s\right)$$

$T_s$ is 300 everywhere, so the diffusion term contributes nothing. $T_s = T_e = T_p = 300$ in the nonmetal domain, so the sources contribute nothing. Therefore, $T_s$ never changes.

Several issues beyond that:

• You set both C_s and K_s to zero in the gold and chromium domains. Similarly with C_e, C_p, and K_e in the nonmetal domain. This gives undefined behavior. You're basically saying that the diffusion coefficient is 0/0 in that part of the simulation. Some solvers (like SciPy PCG) tolerate this and return an answer (although it may be nonsense). Other solvers (SciPy LU, PETSc PCG) raise a singular matrix error. I recommend you set all of the C_i to something nonzero in all domains.
• By defining all of your domains using <= and >=, it is ambiguous which cells belong to which domain. Better to say, e.g.,

  gold = (x <= Lx / 4.0)
  chromium = (x > Lx / 4.0 ) & ( x<= Lx/2.0)
  metal = (x<= Lx/2.0)
  nonmetal = (x> Lx/2.0)

  Otherwise, for example, a cell with center exactly at Lx / 4.0 will be defined as both gold and chromium.
• Your "Interface Resistance Condition(s)" aren't represented in your code. Where in the domain is $L$?

[Sandip433]

Sir, (i) the source terms are $\sqrt{\frac{4 \log(2)}{\pi}} \left(\frac{{(1-R)I_0}}{{t_{p}\delta}}\right)\exp\left(-\frac{x}{\delta}-4\log(2)\left(\frac{{\text{time}-2t_{p}}}{{t_{p}}}\right)^2\right)$ which is mathematical expression of Laser energy irradiation. The Laser source gives energy to metal electrons and electrons give energy to phonon and at the bottom of metal layer , both electrons and phonon give energy to non-metal substrate which has only single temperature characteristic.
So the source of energy for non-metal substrate are electrons and phonon.
(ii) I set both C_s and K_s to zero in the gold and chromium domains because the variable does not exist in metal region therefore C_s and K_s do not exist in metal region . Similarly I set C_e, C_p, and K_e have value 0 in the nonmetal domain. But according to your suggestion I changed that.
(ii) I also corrected the domains.
(iv) I wrote the equation of non-metal region
$\frac{\partial C_s T_s}{\partial t} = \nabla\cdot\left(k_s\nabla T_s\right) + M \cdot {\tt nonmetal} \left(T_e - T_s + T_p - T_s\right)$
because at the end of metal layer , the outward flux in terms of electrons and phonon temperature gradients multiplied by respective thermal conductivity is equal to the difference of temperatures between electrons and non-metal ( and phonon and non-metal) multiplied by resistance term(constant). Since writing $\cdot\left(k_e\nabla T_e\right) = S$ (source term) in first equation works , so I wrote the eqns for non-metal regions according to that, replacing flux terms with temperature difference. I don't know better idea to represent Interface Resistance Condition. Interface Resistance Condition means the non-equilibrium in temperature between metal and non-metal region i.e. at x = Lx/2
the updated code still does not do better

from fipy import *

from fipy.solvers import LinearPCGSolver
Lx=2.0e-7 
dx=1.0e-10
nx=int(Lx/dx) 
T=100e-12
dt=4e-15
I_0 = 1100 
R = 0.982
delta = 12.23e-9
t_p = 100e-15
m=Grid1D(nx=nx, dx=dx)
x = m.cellCenters[0]
time = Variable(0.)
gold = (x <= Lx / 4.0)
chromium = (x > Lx / 4.0 ) & ( x<= Lx/2.0)
metal = (x<= Lx/2.0)
nonmetal = (x> Lx/2.0)
S = numerix.sqrt(4*numerix.log(2)/numerix.pi)*((1-R)*I_0/(t_p*delta))*numerix.exp(-(x/delta)-4*numerix.log(2)*((time-2*t_p)/t_p)**2)
T_e=CellVariable(name="Electron Temperature", mesh=m, value=300.0, hasOld=True )
T_p=CellVariable(name="Phonon Temperature", mesh=m, value=300.0, hasOld=True )
T_s=CellVariable(name="Non-metal Temperature", mesh=m, value=300.0, hasOld=True )
C_e = (68 * gold + 194 * chromium + 1 * nonmetal) * T_e
C_p = (2.5e6 * gold + 3.3e6 * chromium + 1 * nonmetal)
C_s = 2e6 
G_g = 0.21e17*((1.18e7/1.25e11)*(T_e+T_p)+1)
G_c = 4.2e17*((7.9e7/13.4e11)*(T_e+T_p)+1)
K_eg = 318*(1.25e11*T_e/(1.18e7*T_e**2+1.25e11*T_p))
K_ec = 95*(13.4e11*T_e/(7.9e7*T_e**2+13.4e11*T_p))
K_e = K_eg * gold + K_ec * chromium + 1 * nonmetal
K_s = 10 
G = G_g * gold + G_c * chromium + 1 * nonmetal
M=1e10 # taking R_es=R_ps=1e10
eq_e=TransientTerm ( coeff= C_e , var=T_e) == DiffusionTerm (coeff=K_e, var=T_e) - ImplicitSourceTerm (coeff=G, var=T_e) + ImplicitSourceTerm (coeff=G, var=T_p) + S * gold
eq_p=TransientTerm ( coeff= C_p , var=T_p) == DiffusionTerm (coeff=(0.01 * K_e ) , var=T_p) + ImplicitSourceTerm (coeff=G, var=T_e) - ImplicitSourceTerm (coeff=G, var=T_p)
eq_s=TransientTerm ( coeff= C_s , var=T_s) == DiffusionTerm (coeff=K_s, var=T_s) + ImplicitSourceTerm (coeff=  M * nonmetal, var=T_e) - ImplicitSourceTerm (coeff= M * nonmetal, var=T_s) + ImplicitSourceTerm (coeff= M * nonmetal, var=T_p) - ImplicitSourceTerm (coeff=M * nonmetal, var=T_s) 
eq=eq_e & eq_p & eq_s
solver = LinearPCGSolver(precon=None, tolerance=1e-10, iterations=2000)
while time()<T:  
    T_e.updateOld()
    T_p.updateOld()
    T_s.updateOld()
    for i in range(10):
         res=eq.sweep(dt=dt, solver=solver)
    time.setValue(time() + dt)

[guyer]

Think about what can make $T_s$ change:

• $S$ is only applied in the gold phase and is only a source of $T_e$, so it can't make $T_s$ change.
• $T_s$ is initially uniform, so there is no diffusion of this field.
• You say that $T_s$ doesn't exist in the gold or chromium layers and that $T_e$ and $T_p$ don't exist in the nonmetal layer, so there is never a temperature difference that would cause $M \cdot {\tt nonmetal} (T_e - T_s)$ or $M \cdot {\tt nonmetal} (T_p - T_s)$ to give a change in $T_s$

If nothing can cause $T_s$ to change, then it won't change. For electrons or phonons to induce changes in substrate temperature, these fictitious temperatures must exist in the same place somewhere.

I neglected one thing, however. At the face that separates chromium from nonmetal, the diffusion coefficient for $T_e$ and for $T_p$ is the arithmetic mean between the value on the chromium side ($K_{ec}$ or $0.01 K_{ec}$) and the value on the nonmetal side (1, previously 0). This will lead to a diffusive flux of $T_e$ and $T_p$ into the nonmetal cell closest to the chromium layer. With the previous case of $K_e = 0$, these will not diffuse any further into the nonmetal layer, but this slight leakage of $T_e$ and $T_p$ into the chromium layer will enable the $M \cdot {\tt nonmetal} (T_e - T_s + T_p - T_s)$ source of $T_s$ in that first cell, which can then start to diffuse deeper into the nonmetal layer. If you enable a Viewer for $T_s$, you'll see that it does change, although by only a tiny amount. The change is so small because $M/C_s$ is nine orders of magnitude smaller than $G/C_e$.

note: you do not need to set both $C_i$ and $K_i$ to one in the phases where the corresponding temperature "doesn't exist" to keep the equations from being singular. I recommend setting each $C_i$ to a nonzero value everywhere and letting $K_i=0$ in any phase where you don't want the corresponding temperature to conduct. If it doesn't make any sense for electrons and phonons to induce each other in the substrate and they should only induce substratons, then let $G=0$ in the nonmetal layer, too.

[Sandip433]

Yes Sir, the non-metal temperature T_s does change but it is significantly small , it is in 1e-7 range but it should be in range of 1e1 to 1e2 . So I think I wrote third equation wrong. How to implement this conditions
$-\kappa_e \frac{\partial T_e}{\partial x} \bigg| _ {x=L_x}$ = $\frac{1}{R{es}} \bigg( T_e - T_s \bigg) \bigg| _ {x=L_x}$
$-\kappa_p \frac{\partial T_p}{\partial x} \bigg| _ {x=L_x}$ = $\frac{1}{R_{ps}} \bigg( T_p - T_s \bigg) \bigg| _ {x=L_x}$
$-\kappa_s \frac{\partial T_s}{\partial x} \bigg| _ {x=L_x} = \frac{1}{R_{ps}} \bigg( T_p - T_s \bigg) \bigg| _ {x=L_x} + \frac{1}{R{ps}} \bigg( T_p - T_s \bigg) \bigg| _ {x=L_x}$
in this equation $\frac{\partial C_s T_s}{\partial t} = \nabla\cdot\left(k_s\nabla T_s\right)$
instead of $\frac{\partial C_s T_s}{\partial t} = \nabla\cdot\left(k_s\nabla T_s\right) + \frac{1}{R_{es}} \cdot {\tt nonmetal} \left(T_e - T_s \right) + \frac{1}{R_{ps}} \cdot {\tt nonmetal} \left(T_p - T_s \right)$ ?
Here is my updated code
`from fipy import *

from fipy.solvers import LinearPCGSolver
Lx=2.0e-7
dx=1.0e-10
nx=int(Lx/dx)
T=50e-12
dt=4e-15
I_0 = 50
R = 0.359
delta = 115e-9
t_p = 200e-15
m=Grid1D(nx=nx, dx=dx)
x = m.cellCenters[0]
time = Variable(0.)
gold = (x <= 9 Lx / 20.0)
nickel = (x> 9 Lx / 20.0) & ( x <= Lx/2.0)
nonmetal = (x> Lx /2.0 )
S = numerix.sqrt(4 * numerix.log(2) / numerix.pi) * ((1 - R) * I_0 / (t_p * delta * (1 - numerix.exp(-Lx / (2.0 * delta))))) * numerix.exp(-(x / delta) - 4 * numerix.log(2) * ((time - 2 * t_p) / t_p)**2)
T_e=CellVariable(name="Electron Temperature", mesh=m, value=300.0, hasOld=True )
T_p=CellVariable(name="Phonon Temperature", mesh=m, value=300.0, hasOld=True )
T_s=CellVariable(name="Non-metal Temperature", mesh=m, value=300.0, hasOld=True )
C_e = (68 * gold + 1065 nickel + 0.001 * nonmetal) * T_e
C_p = (2.5e6 * gold + 4.1e6 * nickel + 0.001 * nonmetal)
C_s = 2.21e6 #(0.001 * gold + 0.001 * chromium + 2e6 * nonmetal)
G_g = 1.5e16((1.2e7/3.6e11)(T_e+T_p)+1)
G_n = 3.6e17((0.59e7/1.4e11)(T_e+T_p)+1)
Theta_e=T_e/6.42e4
Theta_p=T_p/6.42e4
K_eg = 353((Theta_e*(Theta_e2 + 0.44)*(Theta_e2 + 0.16)1.25)/((Theta_e2 + 0.16Theta_p)(Theta_e**2 + 0.092)0.5))
K_en = 901.4e11T_e/(0.59e7*T_e2 + 1.4e11T_p)
K_e = K_eg * gold + K_en nickel + 0 * nonmetal
K_s = (0 * gold + 0 * nickel + 9.724 * nonmetal)
G = G_g * gold + G_n * nickel + 0 * nonmetal
R_es=1/4.3e-9
R_ps=1/2.06e-8
eq_e=TransientTerm ( coeff= C_e , var=T_e) == DiffusionTerm (coeff=K_e, var=T_e) - ImplicitSourceTerm (coeff=G, var=T_e) + ImplicitSourceTerm (coeff=G, var=T_p) + S * gold
eq_p=TransientTerm ( coeff= C_p , var=T_p) == DiffusionTerm (coeff=(0.01 * K_e ) , var=T_p) + ImplicitSourceTerm (coeff=G, var=T_e) - ImplicitSourceTerm (coeff=G, var=T_p)
eq_s=TransientTerm ( coeff= C_s , var=T_s) == DiffusionTerm (coeff=K_s, var=T_s) + ImplicitSourceTerm (coeff= R_es * nonmetal, var=T_e) - ImplicitSourceTerm (coeff= R_es * nonmetal, var=T_s) + ImplicitSourceTerm (coeff= R_ps * nonmetal, var=T_p) - ImplicitSourceTerm (coeff=R_ps * nonmetal, var=T_s)
eq=eq_e & eq_p & eq_s
results_T_e = []
results_T_p = []
results_T_s = []
t_list=[]
solver = LinearPCGSolver(precon=None, tolerance=1e-10, iterations=2000)
while time()<T:
T_e.updateOld()
T_p.updateOld()
T_s.updateOld()
t_list.append(float(time()))
results_T_e.append(float(T_e[int(Lx/(2.0dx))-1]))
results_T_p.append(float(T_p[int(Lx/(2.0dx))-1]))
results_T_s.append(float(T_s[int(Lx/(2.0*dx))]))
for i in range(1):
res=eq.sweep(dt=dt, solver=solver)
time.setValue(time() + dt)`

[guyer]

Please see the documentation on applying Robin boundary conditions

[Sandip433]

Yes , according to your suggestion I rewrote the code and I understand that the boundary condition is internal robin condition. I rewrote the the boundary condition as such $\frac{1}{R_{es}} \bigg( T_s\bigg) \bigg| _ {x=L_x} + \frac{1}{R_{ps}} \bigg( T_s \bigg) \bigg| _ {x=L_x} -\kappa_s \frac{\partial T_s}{\partial x} \bigg| _ {x=L_x} = \frac{1}{R_{es}} \bigg( T_e \bigg) \bigg| _ {x=L_x} + \frac{1}{R{ps}} \bigg( T_p \bigg) \bigg| _ {x=L_x}$

from fipy import *

from fipy.solvers import LinearPCGSolver
Lx=2.0e-7 
dx=1.0e-10
nx=int(Lx/dx) 
T=30e-12
dt=4e-15
I_0 = 50 
R = 0.359
delta = 115e-9
t_p = 200e-15
m=Grid1D(nx=nx, dx=dx)
x = m.cellCenters[0]
X = m.faceCenters[0]
time = Variable(0.)
gold = (x < Lx / 2.0)
interface = ( X >= Lx/2.0)
nonmetal = (x > Lx /2.0 )
S = numerix.sqrt(4*numerix.log(2)/numerix.pi)*((1-R)*I_0/(t_p*delta*(1-numerix.exp(-Lx/(2.0*delta)))))*numerix.exp(-(x/delta)-4*numerix.log(2)*((time-2.0*t_p)/t_p)**2)
T_e=CellVariable(name="Electron Temperature", mesh=m, value=300.0, hasOld=True )
T_p=CellVariable(name="Phonon Temperature", mesh=m, value=300.0, hasOld=True )
T_s=CellVariable(name="Non-metal Temperature", mesh=m, value=300.0, hasOld=True )
C_e = 68 * T_e
C_p = 2.5e6 
C_s = 1.6e6 
G_g = 1.5e16*((1.2e7/3.6e11)*(T_e+T_p)+1)
Theta_e=T_e/6.42e4
Theta_p=T_p/6.42e4
K_eg = 353*((Theta_e*(Theta_e**2 + 0.44)*(Theta_e**2 + 0.16)**1.25)/((Theta_e**2 + 0.16*Theta_p)*(Theta_e**2 + 0.092)**0.5))
K_e = K_eg * gold   + 0 * nonmetal
K_s = (0 * gold   + 148 * nonmetal)
G = G_g * gold   + 0 * nonmetal
R_es=1/3e-9
R_ps=1e9
a_s = FaceVariable(mesh=m, value=0.0, rank=1)
a_s.setValue((R_es + R_ps)*m.faceNormals, where=interface)
b_s = FaceVariable(mesh=m, value=0.0, rank=0)
b_s.setValue(K_s.faceValue, where=interface)
g_s = FaceVariable(mesh=m, value=0.0, rank=0)
g_s.setValue(R_es*T_e.faceValue + R_ps*T_p.faceValue, where=interface)
eq_e=TransientTerm ( coeff= C_e , var=T_e) == DiffusionTerm (coeff=K_e, var=T_e) - ImplicitSourceTerm (coeff=G, var=T_e) + ImplicitSourceTerm (coeff=G, var=T_p) + S * gold
eq_p=TransientTerm ( coeff= C_p , var=T_p) == DiffusionTerm (coeff=(0.01 * K_e ) , var=T_p) + ImplicitSourceTerm (coeff=G, var=T_e) - ImplicitSourceTerm (coeff=G, var=T_p) 
eq_s=TransientTerm ( coeff= C_s , var=T_s) == DiffusionTerm (coeff=K_s, var=T_s) + PowerLawConvectionTerm(coeff=a_s, var=T_s) + DiffusionTerm(coeff=b_s, var=T_s) + (g_s * interface * m.faceNormals).divergence 
eq=eq_e & eq_p & eq_s
results_T_e = []
results_T_p = []
results_T_s = []
t_list=[]                                               
solver = LinearPCGSolver(precon=None, tolerance=1e-10, iterations=2000)
progress_bar = tqdm(total=int(T/dt), desc='Solving PDE')
while time()<T:  
    T_e.updateOld()
    T_p.updateOld()
    T_s.updateOld()
    t_list.append(float(time()))
    results_T_e.append(float(T_e.faceValue[int(Lx/(2.0*dx))])) 
    results_T_p.append(float(T_p.faceValue[int(Lx/(2.0*dx))])) 
    results_T_s.append(float(T_s.faceValue[int(Lx/(2.0*dx))]))
    for i in range(10):
        res=eq.sweep(dt=dt, solver=solver)
    time.setValue(time() + dt)

and it works fine but when I divide the domain into three parts 90nm gold, 10nm nickel and 100nm nonmetal then the temperature of nonmetal is too much higher than electron temperature

from fipy import *
from fipy.solvers import LinearPCGSolver
Lx=2.0e-7 
dx=1.0e-10
nx=int(Lx/dx) 
T=100e-12
dt=4e-15
I_0 = 50 
R = 0.359
delta = 115e-9
t_p = 200e-15
m=Grid1D(nx=nx, dx=dx)
x = m.cellCenters[0]
X = m.faceCenters[0]
time = Variable(0.)
gold = (x < 9*Lx / 20.0)
nickel =  (x >= 9*Lx / 20.0) & ( x <= Lx / 2.0)
interface = ( X >= Lx/2.0)
nonmetal = (x > Lx /2.0 )
S = numerix.sqrt(4*numerix.log(2)/numerix.pi)*((1-R)*I_0/(t_p*delta*(1-numerix.exp(-Lx/(2.0*delta)))))*numerix.exp(-(x/delta)-4*numerix.log(2)*((time-2.0*t_p)/t_p)**2)
T_e=CellVariable(name="Electron Temperature", mesh=m, value=300.0, hasOld=True )
T_p=CellVariable(name="Phonon Temperature", mesh=m, value=300.0, hasOld=True )
T_s=CellVariable(name="Non-metal Temperature", mesh=m, value=300.0, hasOld=True )
C_e = ( 68 * T_e * gold + 1065 * T_e * nickel )  + 1 * nonmetal
C_p = 2.5e6 * gold + 4.1e6 * nickel + 1 * nonmetal 
C_s = 2.21e6 
G_g = 0.21e17*((1.18e7/1.25e11)*(T_e+T_p)+1)
G_n = 3.6e17*((0.59e7/1.4e11)*(T_e+T_p)+1)
K_eg = 318*(1.25e11*T_e/(1.18e7*T_e**2+1.25e11*T_p))
K_en = 90*(1.4e11*T_e/(0.59e7*T_e**2+1.4e11*T_p))
K_e = K_eg * gold + K_en * nickel + 0 * nonmetal
G = G_g * gold + G_n * nickel + 0 * nonmetal  
K_s = 0 * gold  + 0 * nickel  + 9.724 * nonmetal   
R_es=5e9
R_ps=1/6.7e-9
a_s = FaceVariable(mesh=m, value=0.0, rank=1)
a_s.setValue((R_es + R_ps)*m.faceNormals, where=interface)
b_s = FaceVariable(mesh=m, value=0.0, rank=0)
b_s.setValue(K_s.faceValue, where=interface)
g_s = FaceVariable(mesh=m, value=0.0, rank=0)
g_s.setValue(R_es*T_e.faceValue + R_ps*T_p.faceValue, where=interface)
eq_e=TransientTerm ( coeff= C_e , var=T_e) == DiffusionTerm (coeff=K_e, var=T_e) - ImplicitSourceTerm (coeff=G, var=T_e) + ImplicitSourceTerm (coeff=G, var=T_p) + S * gold
eq_p=TransientTerm ( coeff= C_p , var=T_p) == DiffusionTerm (coeff=(0.01 * K_e ) , var=T_p) + ImplicitSourceTerm (coeff=G, var=T_e) - ImplicitSourceTerm (coeff=G, var=T_p) 
eq_s=TransientTerm ( coeff= C_s , var=T_s) == DiffusionTerm (coeff=K_s, var=T_s) + PowerLawConvectionTerm(coeff=a_s, var=T_s) + DiffusionTerm(coeff=b_s, var=T_s) + (g_s * interface * m.faceNormals).divergence 
eq=eq_e & eq_p & eq_s
results_T_e = []
results_T_p = []
results_T_s = []
t_list=[]                                               
solver = LinearPCGSolver(precon=None, tolerance=1e-10, iterations=2000)
progress_bar = tqdm(total=int(T/dt), desc='Solving PDE')
while time()<T:  
    T_e.updateOld()
    T_p.updateOld()
    T_s.updateOld()
    t_list.append(float(time()))
    results_T_e.append(float(T_e.faceValue[int(Lx/(2.0*dx))])) 
    results_T_p.append(float(T_p.faceValue[int(Lx/(2.0*dx))])) 
    results_T_s.append(float(T_s.faceValue[int(Lx/(2.0*dx))]))
    for i in range(10):
        res=eq.sweep(dt=dt, solver=solver)
    time.setValue(time() + dt)

Did I write the robin condition correctly? If yes then where is the problem?

[guyer]

Your definition of g_s freezes the values of T_e and T_p from the beginning. I would write

R_es_s = FaceVariable(mesh=m, value=0.0, rank=1)
R_es_s.setValue(R_es*m.faceNormals, where=interface)
R_ps_s = FaceVariable(mesh=m, value=0.0, rank=1)
R_ps_s.setValue(R_es*m.faceNormals, where=interface)
  :
  :
eq_s=(TransientTerm ( coeff= C_s , var=T_s)
      == DiffusionTerm (coeff=K_s, var=T_s)
      + PowerLawConvectionTerm(coeff=a_s, var=T_s)
      + DiffusionTerm(coeff=b_s, var=T_s)
      - PowerLawConvectionTerm(coeff=R_es_s, var=T_e)
      - PowerLawConvectionTerm(coeff=R_ps_s, var=T_p))

[Sandip433]

Other than that I have some questions regarding fipy use.
(i) When I am implementing internal condition I can do that with the help of mask idea but in fipy documentation the concept Largevalue is shown. When should I use Largevalue? I applied internal Robin Condition without Largevalue concept . Is it proper way to do that ?
(ii) Instead of this internal boundary condition
$\left(-k_s\nabla T_s\right)\bigg| _ {x=L_x/2} = \frac{1}{R_{es}} \bigg( T_e - T_s \bigg) \bigg| _ {x=L_x/2} + \frac{1}{R{ps}} \bigg( T_p - T_s \bigg) \bigg| _ {x=L_x/2}$
I want to use
$\nabla\cdot\left(k_s\nabla T_s\right)\bigg| _ {x=L_x/2} = \frac{1}{R_{es}} \bigg( T_e - T_s \bigg) \bigg| _ {x=L_x/2} + \frac{1}{R{ps}} \bigg( T_p - T_s \bigg) \bigg| _ {x=L_x/2}$
First condition can be applied through Robin Condition . I tried to apply second condition but T_s did not change .
DiffusionTerm(coeff=largeValue * mask) - ImplicitSourceTerm(mask * largeValue * ConditionValue )
There is two question regarding this.
(a) diffusion term can take facevariable coefficient but source term cannot . So I have to define mask with facevariable in diffuson coefficient and with cellvariable in source term but in a single equation fipy does not allow that.
(b) In case of Robin condition, the whole condition is applied after taking gradient or divergent , how to do that to diffusion term?
(iii) When I am solving Two Temperature Model(two coupled non-linear partial differential equations) , the equations give stable result with only 1 sweep. But if I increase the sweep number to 2 the basic result is stable but not in same range like 1 sweep (temperature for 1 sweep is 2200K and 2 sweep is 1400K). After 20 or 21 sweep the result becomes stable and gives almost identical result(~1600K). When I compare the results with other results solved with FDM or other schemes in some journal paper I see the resemblance of those results with my results I got with 1 sweep, not with 20/21 sweep. This confuses me.

[guyer]

(i) Fundamentally, the cell-centered finite volume discretization that FiPy uses is about solving the divergence equation; values in cells change because of divergences of fluxes into and out of the cell or due to sources and sinks within the cell. This gets broken down into a linear algebra matrix that relates the values of cells and their neighbors. You use a largeValue to overwhelm the matrix solution for the values in the cells of interest because you don't care about their neighbors; the cell should have a specific value or the gradient between cells should have a specific value. A Robin condition is not about the value in the cell but about the flux at the face between cells. A flux may be related to the gradient, but it is not the gradient. In the Robin condition, you are replacing the flux used in the DiffusionTerm or ConvectionTerm with the flux for the Robin condition and you don't need largeValue.

[guyer]

(ii) In the expression
$$\nabla\cdot\left(k_s\nabla T_s\right)\bigg| _ {x=L_x/2} = \frac{1}{R_{es}} \bigg( T_e - T_s \bigg) \bigg| _ {x=L_x/2} + \frac{1}{R{ps}} \bigg( T_p - T_s \bigg) \bigg| _ {x=L_x/2}$$
the flux $k_s\nabla T_s$ occurs at a face bounding a cell, but the divergence $\nabla\cdot\left(k_s\nabla T_s\right)$ occurs at the cell. $\frac{1}{R_{es}} \bigg( T_e - T_s \bigg)$ and $\frac{1}{R{ps}} \bigg( T_p - T_s \bigg)$ also occur at a cell, so there's no inconsistency. This isn't a boundary condition, though, because cells aren't boundaries.

If you accept that this condition applies over a cell volume, then you want to replace the DiffusionTerm in the appropriate cells with this other condition, e.g.,

mask = cells where L_x/2 # this doesn't make sense. L_x/2 is a point, but cells occupy a volume

eq_s = (TransientTerm ( coeff= C_s , var=T_s)
        == DiffusionTerm (coeff=K_s * ~mask, var=T_s)  # turns off K_s in cells of interest
        + ImplicitSourceTerm(mask * R_es, var=T_e) # turns on temperature jumps in cells of interest
        - ImplicitSourceTerm(mask * R_es, var=T_s)
        + ImplicitSourceTerm(mask * R_ps, var=T_p)
        - ImplicitSourceTerm(mask * R_ps, var=T_s))

[guyer]

(iii) This doesn't seem to be related to the other questions. Please post a separate question and include the code you're talking about.

[Sandip433]

Sir, after I tried many variations this is the code I wrote for
Interface.Resistance.pdf

from fipy import *
from fipy.solvers import LinearPCGSolver
import matplotlib.pyplot as plt
import scienceplots
plt.style.use(['science','no-latex','notebook'])
from tqdm import tqdm
Lx=2.0e-7 
dx=1.0e-10
nx=int(Lx/dx) 
T=30e-12
dt=4e-15
I_0 = 50 
R = 0.359
delta = 115e-9
t_p = 200e-15
m=Grid1D(nx=nx, dx=dx)
x = m.cellCenters[0]
X = m.faceCenters[0]
time = Variable(0.)
gold = (x <= 9.0 * Lx/20.0 )
nickel = (x > 9.0 * Lx/20.0 ) & (x <= Lx/2.0 )
nonmetal = (x > Lx/2.0 )
interface = ( X >= Lx/2.0)
S = numerix.sqrt(4*numerix.log(2)/numerix.pi)*((1-R)*I_0/(t_p*delta*(1-numerix.exp(-Lx/(2.0*delta)))))*numerix.exp(-(x/delta)-4*numerix.log(2)*((time-2.0*t_p)/t_p)**2)
T_e=CellVariable(name="Electron Temperature", mesh=m, value=300.0, hasOld=True )
T_p=CellVariable(name="Phonon Temperature", mesh=m, value=300.0, hasOld=True )
T_s=CellVariable(name="Non-metal Temperature", mesh=m, value=300.0, hasOld=True )
C_e = 68 * T_e * gold + 1065 * T_e * nickel  + 1065 * T_e * nonmetal
C_p = 2.5e6 * gold + 4.1e6 * nickel + 4.1e6 * nonmetal
C_s = 2.21e6 
G_g = 1.5e16*((1.2e7/3.6e11)*(T_e+T_p)+1)
G_n = 3.6e17*((0.59e7/1.4e11)*(T_e+T_p)+1)
Theta_e=T_e/6.42e4
Theta_p=T_p/6.42e4
K_eg = 353*((Theta_e*(Theta_e**2 + 0.44)*(Theta_e**2 + 0.16)**1.25)/((Theta_e**2 + 0.16*Theta_p)*(Theta_e**2 + 0.092)**0.5))
K_en = 90*(1.4e11*T_e/(0.59e7*T_e**2+1.4e11*T_p))
K_e = K_eg * gold + K_en * nickel + 0 * nonmetal
G =  G_g * gold + G_n * nickel + 0 * nonmetal  
K_s = 0 * gold + 0 * nickel  + 9.724 * nonmetal   
R_es=1/2e-10
R_ps=1/6.7e-9
k_s = FaceVariable(mesh=m, value=0.0, rank=0)
k_s.setValue(K_s.faceValue, where=interface)
R_es_s = FaceVariable(mesh=m, value=0.0, rank=1)
R_es_s.setValue(R_es*m.faceNormals, where=interface)
R_ps_s = FaceVariable(mesh=m, value=0.0, rank=1)
R_ps_s.setValue(R_es*m.faceNormals, where=interface)
eq_e=TransientTerm ( coeff= C_e , var=T_e) == DiffusionTerm (coeff=K_e, var=T_e) - ImplicitSourceTerm (coeff=G, var=T_e) + ImplicitSourceTerm (coeff=G, var=T_p) + S * gold
eq_p=TransientTerm ( coeff= C_p , var=T_p) == DiffusionTerm (coeff=(0.01 * K_e ) , var=T_p) + ImplicitSourceTerm (coeff=G, var=T_e) - ImplicitSourceTerm (coeff=G, var=T_p) 
eq_s=TransientTerm ( coeff= C_s , var=T_s) == DiffusionTerm (coeff=K_s, var=T_s) + DiffusionTerm (coeff=k_s, var=T_s)  + PowerLawConvectionTerm(coeff=R_es_s, var=T_e) + PowerLawConvectionTerm(coeff=R_ps_s, var=T_p) - PowerLawConvectionTerm(coeff=R_es_s, var=T_s) - PowerLawConvectionTerm(coeff=R_ps_s, var=T_s)
eq=eq_e & eq_p & eq_s
results_T_e = []
results_T_p = []
results_T_s = []
t_list=[]                                               
solver = LinearPCGSolver(precon=None, tolerance=1e-10, iterations=2000)
progress_bar = tqdm(total=int(T/dt), desc='Solving PDE')
while time()<T:  
    T_e.updateOld()
    T_p.updateOld()
    T_s.updateOld()
    t_list.append(float(time()))
    results_T_e.append(float(T_e.faceValue[int(Lx/(2.0*dx) + 1)])) 
    results_T_p.append(float(T_p.faceValue[int(Lx/(2.0*dx) + 1)])) 
    results_T_s.append(float(T_s.faceValue[int(Lx/(2.0*dx) + 1)]))
    for i in range(10):
        res=eq.sweep(dt=dt, solver=solver)
    time.setValue(time() + dt)
    progress_bar.update()
progress_bar.close()
TSVViewer(vars=( T_e.faceValue,T_p.faceValue,T_s.faceValue)).plot(filename="123 2.txt")
plt.plot(t_list,results_T_e)
plt.plot(t_list,results_T_p)
plt.plot(t_list,results_T_s)
plt.xlabel('X (time)')
plt.ylabel('Y (Temperature)')
plt.title('Temperature')
plt.savefig('Interface 2.png', dpi=1200)
plt.show()

But there are still some problems and questions remained.
(i) in this code the temperature of substrate changes only in interface site (X==Lx/2.0) but not in subsequent region diffusively
(ii) To apply flux coming from electron and phonon do I have to apply robin condition in each equation or just in substrate equation the sum of two fluxes coming from electron and phonon? If I apply in electron and phonon equation the equations blow up.
(iii)when we apply robin condition why we apply under divergence with diffusion and convection term ? why not directly with convection term and source term?
(iv)when we apply robin condition , by the definition of b in the fipy documentation the diffusion term is doubled. Is that right?

b = FaceVariable(mesh=mesh, value=K_s, rank=0)
b.setValue(0., where=mask)
g = FaceVariable(mesh=mesh, value=..., rank=0)
eqn = (TransientTerm() == DiffusionTerm(coeff=K_s) + PowerLawConvectionTerm(coeff=a) + DiffusionTerm(coeff=b) + (g * mask * mesh.faceNormals).divergence)

(v)when I mention the interface the proper region is X==Lx/2.0 but this condition does not give any result but if I mention X>=Lx/2.0 then some changes happen. why?
(vi) problems regarding sign.
when we apply robin condition $\cdot\left(a.T_e + k_e\nabla T_e\right) = g$ we write
TransientTerm() == PowerLawConvectionTerm(coeff=a) + DiffusionTerm(coeff=K_e) + (g * mask * mesh.faceNormals).divergence)
but when we write internal fixed gradient we write
(TransientTerm() == DiffusionTerm(coeff=Gamma)+ DiffusionTerm(coeff=largeValue * mask) - ImplicitSourceTerm(mask * largeValue * gradient * mesh.faceNormals).divergence)
why the signs do not matter?
(vii)in fipy documentation it said that The parameter largeValue must be chosen to be large enough to completely dominate the matrix diagonal and the RHS vector in cells that are masked but how to determine that this value is large enough? Do we have to experiment or is there any proper rule?
(viii)in fipy documentation , it said that An expression like the heat flux convection boundary condition −𝑘∇𝑇 · ˆ𝑛 = ℎ(𝑇 −𝑇∞) can be put in the form of the Robin condition used above by letting ⃗𝑎 ≡ ℎˆ𝑛, 𝑏 ≡ 𝑘, and 𝑔 ≡ ℎ𝑇∞. but instead of this condtion if we have to mention −𝑘∇𝑇 · ˆ𝑛=function of variables how to do that?
(ix)if we have to mention fixed flux source −𝑘∇𝑇=constant in mesh facesLeft how to do that , by facegrad.constrain(constant/𝑘,mesh.facesLeft)
or (constant * mesh.facesLeft * mesh.faceNormals).divergence), constant=facevariable rank 1

[guyer]

But there are still some problems and questions remained.

(i) in this code the temperature of substrate changes only in interface site (X==Lx/2.0) but not in subsequent region diffusively

That's not true. After just a few time steps, this is what your substrate temperature looks like:

(ii) To apply flux coming from electron and phonon do I have to apply robin condition in each equation or just in substrate equation the sum of two fluxes coming from electron and phonon? If I apply in electron and phonon equation the equations blow up.

You need to apply the boundary conditions to each field where they apply. Based on your PDF, you have a Robin condition for each of $T_e$, $T_p$, and $T_s$, so apply them to each.

(iii)when we apply robin condition why we apply under divergence with diffusion and convection term ? why not directly with convection term and source term?

I don't understand what you're asking. Please show math and/or code to illustrate what you're talking about.

(iv)when we apply robin condition , by the definition of b in the fipy documentation the diffusion term is doubled. Is that right?

b = FaceVariable(mesh=mesh, value=K_s, rank=0)
b.setValue(0., where=mask)
g = FaceVariable(mesh=mesh, value=..., rank=0)
eqn = (TransientTerm() == DiffusionTerm(coeff=K_s) + PowerLawConvectionTerm(coeff=a) + DiffusionTerm(coeff=b) + (g * mask * mesh.faceNormals).divergence)

No, because you shouldn't have DiffusionTerm(coeff=K_s) in your equation.

(v)when I mention the interface the proper region is X==Lx/2.0 but this condition does not give any result but if I mention X>=Lx/2.0 then some changes happen. why?

X==Lx/2.0 isn't a region, it's a location. If that location does not correspond exactly to the location of a face center, then there will be no faces described by that location. Try

interface = ( X > (Lx - dx)/2.0) & ( X < (Lx + dx)/2.0)

(vi) problems regarding sign. when we apply robin condition ⋅(a.Te+ke∇Te)=g we write TransientTerm() == PowerLawConvectionTerm(coeff=a) + DiffusionTerm(coeff=K_e) + (g * mask * mesh.faceNormals).divergence) but when we write internal fixed gradient we write (TransientTerm() == DiffusionTerm(coeff=Gamma)+ DiffusionTerm(coeff=largeValue * mask) - ImplicitSourceTerm(mask * largeValue * gradient * mesh.faceNormals).divergence) why the signs do not matter? tial

What do you mean the signs don't matter? Those aren't doing the same thing. For an internal fixed gradient

$$ \begin{align} \frac{\partial\phi}{\partial t} &= \nabla\cdot\left(\Gamma_0 \nabla\phi\right) \\ \left.\nabla\phi\right|_\mathtt{mask} &= \mathtt{gradient} \end{align} $$

becomes

$$ \begin{align} \frac{\partial\phi}{\partial t} &= \nabla\cdot\left(\Gamma_0 \cdot\mathtt{\sim mask}\ \nabla\phi\right) {} + \nabla\cdot\left(\mathtt{largeNumber}\cdot\mathtt{mask}\ \nabla\phi\right) {} - \nabla\cdot\left(\mathtt{largeNumber}\cdot\mathtt{mask} \cdot \mathtt{gradient}\right) \end{align} $$

When mask is false, this reduces to

$$ \begin{align} \frac{\partial\phi}{\partial t} &= \nabla\cdot\left(\Gamma_0 \nabla\phi\right) \end{align} $$

When mask is true, this reduces to

$$ \begin{align} 0 &= \nabla\cdot\left(\mathtt{largeNumber}\ \nabla\phi\right) {} - \nabla\cdot\left(\mathtt{largeNumber}\cdot \mathtt{gradient}\right) \\ \nabla\cdot\left(\mathtt{largeNumber}\ \nabla\phi\right) &= \nabla\cdot\left(\mathtt{largeNumber}\cdot \mathtt{gradient}\right) \\ \nabla\phi &= \mathtt{gradient} \end{align} $$

(vii)in fipy documentation it said that The parameter largeValue must be chosen to be large enough to completely dominate the matrix diagonal and the RHS vector in cells that are masked but how to determine that this value is large enough? Do we have to experiment or is there any proper rule?

Experiment but, basically, largeNumber needs to be much bigger than your diffusion coefficient.

(viii)in fipy documentation , it said that An expression like the heat flux convection boundary condition −𝑘∇𝑇 · ˆ𝑛 = ℎ(𝑇 −𝑇∞) can be put in the form of the Robin condition used above by letting ⃗𝑎 ≡ ℎˆ𝑛, 𝑏 ≡ 𝑘, and 𝑔 ≡ ℎ𝑇∞. but instead of this condtion if we have to mention −𝑘∇𝑇 · ˆ𝑛=function of variables how to do that?

It depends on the function. If it's a linear function of T, then subsitute appropriately for $\hat{n}\cdot\vec{a} T$. If it's not, then substitute for $g$. These are just algebraic manipulations.

(ix)if we have to mention fixed flux source −𝑘∇𝑇=constant in mesh facesLeft how to do that , by facegrad.constrain(constant/𝑘,mesh.facesLeft) or (constant * mesh.facesLeft * mesh.faceNormals).divergence), constant=facevariable rank 1

Either should work, but the latter is more physical. Treating fluxes and gradients like they're equivalent is the source of everything wrong in the world.

In the future, can I kindly request that you ask questions as they come up? It's a lot of work to parse through a wall of nine questions after two weeks.

[Sandip433]

I am very sorry for asking questions after long time. But I always want to use your suggestion to make my code work and experiment first for some time.

That's not true. After just a few time steps, this is what your substrate temperature looks like.

I meant that change in interface layer is higher than subsequent facecenter. The change is not consistant as it should be.

You need to apply the boundary conditions to each field where they apply. Based on your PDF, you have a Robin condition for each of T_e,T_p and T_s so apply them to each.

If I apply the robin condition to equation of T_e and T_p the solutions do not converge at all.

I don't understand what you're asking. Please show math and/or code to illustrate what you're talking about.

I meant that when we apply robin condition we take divergence in both sides of $\left(a.T_e + k_e\nabla T_e\right) = g$ to $\nabla\left (a.T_e + k_e\nabla T_e\right) =\nabla g$ and then $\nabla\left (a.T_e) + \nabla (k_e\nabla T_e\right) =\nabla g$ but why not in covection term and source term to represent robin condition?

No, because you shouldn't have DiffusionTerm(coeff=K_s) in your equation

But the main equation contains transient term and diffusion term. Then robin condition comes . Why do we need to exclude diffusion term of main equation?

It depends on the function. If it's a linear function of T, then subsitute appropriately for n.aT. If it's not, then substitute for g. These are just algebraic manipulations.

I want to simulate these equations.
Interface.Resistance.pdf
So the robin condition should be in this form $\left(-k_s\nabla T_s\right)\bigg| _ {x=L_x/2} = \frac{1}{R_{es}} \bigg( T_e - T_s \bigg) \bigg| _ {x=L_x/2} + \frac{1}{R{ps}} \bigg( T_p - T_s \bigg) \bigg| _ {x=L_x/2}$
or $\frac{1}{R_{es}} \bigg( T_s\bigg) \bigg| _ {x=L_x} + \frac{1}{R_{ps}} \bigg( T_s \bigg) \bigg| _ {x=L_x} -\kappa_s \frac{\partial T_s}{\partial x} \bigg| _ {x=L_x} = \frac{1}{R_{es}} \bigg( T_e \bigg) \bigg| _ {x=L_x} + \frac{1}{R{ps}} \bigg( T_p \bigg) \bigg| _ {x=L_x}$
or in this form $\left(-k_s\nabla T_s\right)\bigg| _ {x=L_x/2} = \left(-k_e\nabla T_e\right)\bigg| _ {x=L_x/2} + \left(-k_p\nabla T_p\right)\bigg| _ {x=L_x/2}$
where $\left(-k_e\nabla T_e\right)\bigg| _ {x=L_x/2} = \frac{1}{R_{es}} \bigg( T_e - T_s \bigg) \bigg| _ {x=L_x/2} $
and $\left(-k_p\nabla T_p\right)\bigg| _ {x=L_x/2} = \frac{1}{R{ps}} \bigg( T_p - T_s \bigg) \bigg| _ {x=L_x/2}$ ?
Because I have seen that none of these forms make the substrate change according to this reference picture of substrate temperature change
1.pdf

Sir, according to the Fourier's Law of heat conduction, heat flux in positive direction is mathematically - thermal conductivity * gradient of temperature. But what I try to simulate is that at the end of metal layer the heat flux from electron of metals to substrate depends linearly of temperature difference of electron and substrate, and the heat flux from phonon of metals to substrate depends linearly of temperature difference of phonon and substrate. I know that you have said many times that flux and gradient are not same. So what should I do to represent this flux condition in fipy? Should it be in thermal conductivity * gradients or in diffuson codition ? This is why I asked this question

Instead of this internal boundary condition
$\left(-k_s\nabla T_s\right)\bigg| _ {x=L_x/2} = \frac{1}{R_{es}} \bigg( T_e - T_s \bigg) \bigg| _ {x=L_x/2} + \frac{1}{R{ps}} \bigg( T_p - T_s \bigg) \bigg| _ {x=L_x/2}$
I want to use
$\nabla\cdot\left(k_s\nabla T_s\right)\bigg| _ {x=L_x/2} = \frac{1}{R_{es}} \bigg( T_e - T_s \bigg) \bigg| _ {x=L_x/2} + \frac{1}{R{ps}} \bigg( T_p - T_s \bigg) \bigg| _ {x=L_x/2}$
First condition can be applied through Robin Condition . I tried to apply second condition but T_s did not change .