How to improve first order boundary condition to second order

[Sandip433]

Hi, I was trying to implement heat resistance between metal (100 nm) and non-metal layer(100 nm) , the conditions were
$-\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}$ ( $x= L_x$ is the position of interface )( Internal Robin Condition )
in this equations
$\frac{\partial C_e T_e}{\partial t} = \nabla\cdot\left(k_e\nabla T_e\right) - G(T_e - T_p) + S$
$\frac{\partial C_p T_p}{\partial t} = \nabla\cdot\left(k_p\nabla T_p\right) + G(T_e - T_p)$
$\frac{\partial C_s T_s}{\partial t} = \nabla\cdot\left(k_s\nabla T_s\right)$
First I set up a domain of 200 nm with three variables $T_e, T_p, T_s$ , but I could not apply the internal robin condition to my satisfaction . So I changed the approach to solve the problem. Instead of coupling the equations through heat flux in internal position , I coupled them at boundary position. Here is the code

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=1.0e-7 
dx=1.0e-10
nx=int(Lx/dx) 
T=20e-12
dt=4e-15
I_0 = 50 
R = 0.595
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/10)
nickel = (( x > 9*Lx/10) & ( x < Lx))
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.0 * T_e
C_p = 2.5e6 
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)
G = G_g * gold + G_n * nickel
Theta_e=T_e.faceValue /6.42e4
Theta_p=T_p.faceValue /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.faceValue / (0.59e7 * T_e.faceValue ** 2 + 1.4e11 * T_p.faceValue )
R_es=1/2e-10
R_ps=1/6.7e-9
#K_e = ~m.facesRight * K_eg
#K_s = ~m.facesRight * 9.724
K_e = FaceVariable(mesh=m, value= K_eg * ( X <= 9 * Lx/10) + K_en * (( X > 9 * Lx/10) & ( X < Lx)), rank=0)
K_e.constrain(0.0, where=m.facesRight)
K_s = FaceVariable(mesh=m, value= 9.724, rank=0)
K_s.constrain(0.0, where=m.facesRight)
boundary_interface = (m.facesRight * m.faceNormals).divergence
d_PR = m._cellDistances[m.facesRight.value][0]
A_f = m._faceAreas[m.facesRight.value][0]
boundary_coeff = boundary_interface * A_f / (1.0 + d_PR) 
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 - ImplicitSourceTerm(coeff= R_es * boundary_coeff, var=T_e) + ImplicitSourceTerm(coeff= R_es * boundary_coeff, var=T_s)
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) - ImplicitSourceTerm(coeff= R_ps * boundary_coeff, var=T_p) + ImplicitSourceTerm(coeff= R_ps * boundary_coeff, var=T_s)
eq_s=TransientTerm(coeff=C_s,var=T_s)==DiffusionTerm(coeff=K_s,var=T_s) + ImplicitSourceTerm(coeff= R_es * boundary_coeff, var=T_e) - ImplicitSourceTerm(coeff= R_es * boundary_coeff, var=T_s) + ImplicitSourceTerm(coeff= R_ps * boundary_coeff, var=T_p) - ImplicitSourceTerm(coeff= R_ps * boundary_coeff, 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[int(Lx/dx)-1])) 
    results_T_p.append(float(T_p[int(Lx/dx)-1])) 
    results_T_s.append(float(T_s[int(Lx/dx)-1]))
    for i in range(1):
        res=eq.sweep(dt=dt, solver=solver)
    time.setValue(time() + dt)
    progress_bar.update()
progress_bar.close()
TSVViewer(vars=( T_e,T_p,T_s)).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 sweep 10.png', dpi=1200)
plt.show()
#data = numerix.genfromtxt('123 2.txt', delimiter='\t', skip_header=1)
#data[:, 3] = data[::-1, 3]
#numerix.savetxt('123 2.txt', data, delimiter='\t', fmt='%s')

But previously @wd15 said that this implementation is only first order accurate and he also mentioned that this implementation may have trouble in higher dimension . I want to know that how to correct the boundary condition up to second order?

[guyer]

Related discussions: #1007 #942

But previously @wd15 said that this implementation is only first order accurate and he also mentioned that this implementation may have trouble in higher dimension .

Please provide a reference to that discussion.

[Sandip433]

These conversations helped me to write this code

• https://stackoverflow.com/questions/42769749/general-boundary-conditions
• https://stackoverflow.com/questions/61632193/1d-coupled-transient-diffusion-in-fipy-with-reactive-boundary-condition
• https://stackoverflow.com/questions/46406985/translate-2-d-pdes-and-conditions-into-accurate-fipy-codes