如何在FiPy中添加潜热项？

Question

我一直在尝试使用fipy模拟双温度模型 模型的数学模型：

C_e(∂T_e)/∂t=∇[k_e∇T_e ]-G(T_e-T_ph )+ A(r,t)

C_ph(∂T_ph)/∂t=∇[k_ph∇T_ph] + G(T_e-T_ph)

源应该加热电子{}，然后热通过{}传递到声子{}，当{}达到熔点，例如{}，一些热{}在熔化前作为潜热

这是我的密码：

from fipy.tools import numerix
import  scipy
import fipy
import numpy as np
from fipy import CylindricalGrid1D
from fipy import Variable, CellVariable, TransientTerm, DiffusionTerm, Viewer, LinearLUSolver, LinearPCGSolver, \
    LinearGMRESSolver, ImplicitDiffusionTerm, Grid1D, ImplicitSourceTerm

## Mesh


nr = 50
dr = 1e-7
# r = nr * dr
mesh = CylindricalGrid1D(nr=nr, dr=dr, origin=0)
x = mesh.cellCenters[0]

# Variables
T_e = CellVariable(name="electronTemp", mesh=mesh,hasOld=True)
T_e.setValue(300)
T_ph = CellVariable(name="phononTemp", mesh=mesh, hasOld=True)
T_ph.setValue(300)
G = CellVariable(name="EPC", mesh=mesh)
t = Variable()
# Material parameters
C_e = CellVariable(name="C_e", mesh=mesh)
k_e = CellVariable(name="k_e", mesh=mesh)

C_ph = CellVariable(name="C_ph", mesh=mesh)
k_ph = CellVariable(name="k_ph", mesh=mesh)

C_e = 4.15303 - (4.06897 * numerix.exp(T_e / -85120.8644))
C_ph = 4.10446 - 3.886 * numerix.exp(-T_ph / 373.8)
k_e = 0.1549 * T_e**-0.052
k_ph =1.24 + 16.29 * numerix.exp(-T_ph / 151.57)


G = numerix.exp(21.87 + 10.062 * numerix.log(numerix.log(T_e )- 5.4))

# Boundary conditions
T_e.constrain(300, where=mesh.facesRight)
T_ph.constrain(300, where=mesh.facesRight)

# Source  𝐴(𝑟,𝑡) = 𝑎𝐷(𝑟)𝜏−1 𝑒−𝑡/𝜏   , 𝐷(𝑟) = 𝑆𝑒 exp (−𝑟2/𝜎2)/√2𝜋𝜎2
sig = 1.0e-6
tau = 1e-15
S_e = 35

d_r = (S_e * 1.6e-9 * numerix.exp(-x**2 /sig**2)) / (numerix.sqrt(2. * 3.14 * sig**2))
A_t = numerix.exp(-t/tau)
a = (numerix.sqrt(2. * 3.14)) / (3.14 * sig)

A_r = a * d_r * tau**-1 * A_t



eq0 = (
    TransientTerm(var=T_e, coeff=C_e) == \
    DiffusionTerm(var=T_e, coeff=k_e) - \
    ImplicitSourceTerm(coeff=G, var=T_e) + \
    ImplicitSourceTerm(var=T_ph, coeff=G) + \
    A_r)

eq1 = (TransientTerm(var=T_ph, coeff=C_ph) == DiffusionTerm(var=T_ph, coeff=k_ph) + ImplicitSourceTerm(var=T_e, coeff=G) - ImplicitSourceTerm(coeff=G, var=T_ph))
eq = eq0 & eq1

dt = 1e-18
steps = 7000
elapsed = 0.
vi = Viewer((T_e, T_ph), datamin=0., datamax=2e4)


for step in range(steps):
    T_e.updateOld()
    T_ph.updateOld()
    vi.plot()
    res = 1e100
    dt *= 1.01
    count = 0
    while res > 1:
        res = eq.sweep(dt=dt, underRelaxation=0.5)
        print(t, res)
    t.setValue(t + dt)

据我所知，我可以将潜热作为源项和汇项包含在eq1中，或者在C_ph中添加一个高斯峰，峰中心应该在熔点附近

我不知道哪一个更好、更稳定，我不知道如何实现其中任何一个

请帮帮我

Answer 1

根据评论（请将其编辑成问题），将eq1更改为

eq1 = (TransientTerm(var=T_ph, coeff=C_ph) 
       == DiffusionTerm(var=T_ph, coeff=k_ph)
       + ImplicitSourceTerm(var=T_e, coeff=G) 
       - ImplicitSourceTerm(coeff=G, var=T_ph) 
       + (1/numerix.sqrt(2*numerix.pi * sig2)) * numerix.exp(-(T_ph - 1850)**2 / 2 * sig2)))

它将被显式计算，但只要T_ph更新，它就会更新