Internal Boundary condition of coupled equations

[Sandip433]

I am trying to solve four coupled equations about laser irradiation on bimetallic surface. The equations are the equations are PhysRevB. Equations.pdf I defined the variables `Lx=2.0e-7 #m(depth of consideration)

dx=1.0e-10 #m(grid size) nx=Lx/dx #number of cells in mesh m=Grid1D(nx=nx, dx=dx) x = m.cellCenters[0] time = Variable(0.) T_eg=CellVariable(name="Electron Temperature gold", mesh=m, value=300.0, hasOld=True) T_pg=CellVariable(name="Phonon Temperature gold", mesh=m, value=300.0, hasOld=True) T_ec=CellVariable(name="Electron Temperature chromium", mesh=m, value=300.0, hasOld=True) T_pc=CellVariable(name="Phonon Temperature chromium", mesh=m, value=300.0, hasOld=True) ` But the boundary condition I could not write as required at x= Lx/2

I tried to write first two conditions `T_eg.constrain([T_ec], where = ( x== 1e-7))

T_pg.constrain([T_pc], where = ( x== 1e-7))` Can anybody help?

[guyer]

Please see https://www.ctcms.nist.gov/fipy/documentation/USAGE.html#applying-internal-boundary-conditions

[Sandip433]

After a lot of trial and error, I came up with this code which runs without any error, but I did not get satisfactory result , some variables have undesirable effects. the equations are PhysRevB. Equations.pdf

from fipy import *
import matplotlib.pyplot as plt
import scienceplots
plt.style.use(['science','no-latex','notebook'])
Lx=2.0e-7 #m(depth of consideration)
dx=1.0e-10 #m(grid size)
nx=int(Lx/dx)#number of cells in mesh
I_0 = 1100 
R = 0.982
delta = 12.23e-9
t_p = 100e-15
m=Grid1D(nx=nx, dx=dx) #grid1d does not have top and bottom faces
x = m.cellCenters[0]
time = Variable(0.)

mask = (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_eg=CellVariable(name="Electron Temperature gold", mesh=m, value=300.0, hasOld=True )
T_pg=CellVariable(name="Phonon Temperature gold", mesh=m, value=300.0, hasOld=True )
T_ec=CellVariable(name="Electron Temperature chromium", mesh=m, value=300.0, hasOld=True)
T_pc=CellVariable(name="Phonon Temperature chromium", mesh=m, value=300.0, hasOld=True) 

G_g = 0.21e17 * ((1.18e7/1.25e11) * (T_eg+T_pg)+1)
G_c = 4.2e17 * ((7.9e7/13.4e11) * (T_ec+T_pc)+1)
K_eg = 318 * (1.25e11 * T_eg /(1.18e7 * T_eg ** 2+1.25e11 * T_pg))
K_ec = 95 * (13.4e11 * T_ec /(7.9e7 * T_ec **  2+13.4e11 * T_pc))
T_eg.faceGrad.constrain([0], m.facesLeft)
T_ec.faceGrad.constrain([0], m.facesRight)            
T_pg.faceGrad.constrain([0], m.facesLeft) 
T_pc.faceGrad.constrain([0], m.facesRight) 

eq_e=TransientTerm ( coeff=68*T_eg, var=T_eg) == DiffusionTerm (coeff=K_eg, var=T_eg) - ImplicitSourceTerm (coeff=G_g, var=T_eg) + ImplicitSourceTerm (coeff=G_g, var=T_pg)+S
eq_p=TransientTerm ( coeff=2.5e6, var=T_pg) == DiffusionTerm (coeff=0.01*K_eg, var=T_pg)+ ImplicitSourceTerm (coeff=G_g, var=T_eg) - ImplicitSourceTerm (coeff=G_g, var=T_pg)

eq_a=(TransientTerm (coeff=194 * T_ec, var=T_ec) == DiffusionTerm (coeff=K_ec, var=T_ec) - ImplicitSourceTerm (coeff=G_c, var=T_ec) + ImplicitSourceTerm (coeff=G_c, var=T_pc))- ImplicitSourceTerm (coeff=mask * 1e40, var=T_ec)+ ImplicitSourceTerm (coeff=mask * 1e40  , var=T_eg) - DiffusionTerm (coeff=1e40 * mask * K_ec, var=T_ec) + DiffusionTerm (coeff=1e40 * mask * K_eg, var=T_eg)
eq_t=(TransientTerm (coeff=3.3e6, var=T_pc) == DiffusionTerm (coeff=0.01 * K_ec, var=T_pc) + ImplicitSourceTerm (coeff=G_c, var=T_ec) + ImplicitSourceTerm (coeff=G_c, var=T_pc)) - ImplicitSourceTerm (coeff=mask * 1e40, var=T_pc)+ ImplicitSourceTerm (coeff=mask * 1e40  , var=T_pg) - DiffusionTerm (coeff=1e40 * mask * K_ec, var=T_pc) + DiffusionTerm (coeff=1e40 * mask * K_eg, var=T_pg)
eq=eq_e & eq_p & eq_a & eq_t
results_T_eg = []
results_T_pg = []
results_T_ec = []
results_T_pc = []
t_list=[]
#solver = LinearLUSolver()
while time()<10e-12:
    t_list.append(float(time()))
    results_T_eg.append(float(T_eg.max()))
    results_T_pg.append(float(T_pg.max()))
    results_T_ec.append(float(T_ec.max())) 
    results_T_pc.append(float(T_pc.max()))
    T_eg.updateOld()
    T_pg.updateOld()
    T_ec.updateOld()
    T_pc.updateOld()
    for i in range(1):
        res=eq.sweep(dt=4e-15)
    time.setValue(time() + 4e-15)

plt.plot(t_list,results_T_eg,'red',label="Gold Electron Temperature")
plt.plot(t_list,results_T_pg,'yellow',label="Gold Phonon Temperature")
plt.plot(t_list,results_T_ec,'blue',label="Chromium Electron Temperature")
plt.plot(t_list,results_T_pc,'green',label="Chromium Phonon Temperature")
plt.show()

I tried with + DiffusionTerm(coeff=largeValue * mask) - ImplicitSourceTerm(mask * largeValue * gradient * mesh.faceNormals).divergence) as documented but errors occured. What improvements can be done to this code?

[guyer]

"but errors occured"!?

[guyer]

I don't think you want four equations with four unknowns. You have two equations with two unknowns $T_e$ and $T_p$ and you have two domains with different coefficients for $x < L_1$ and $x > L_1$. If you do that, then the "internal boundary condition" will be satisfied naturally.

[guyer]

E.g.

gold = (x < Lx / 2)
chromium = (x >= Lx / 2)
K_e = K_eg * gold + K_ec * chromium
G = G_g * gold + G_c * chromium

eq_e=TransientTerm ( coeff=(68 * gold + 194 * chromium) * T_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=(2.5e6 * gold + 3.3e6 * chromium), 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 = eq_e & eq_p

[Sandip433]

I am sorry for not clarifying the error, which was cellvariables does not have any attribute 'gradient'. But according to your suggestion the code runs perfectly. Thanks for the help.