如何在 FiPy 中添加潜热项？答案

【问题标题】：How to add latent heat term in FiPy?如何在 FiPy 中添加潜热项？
【发布时间】：2021-09-14 22:27:01
【问题描述】：

我一直在尝试使用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)

源应该加热电子T_e，然后热量通过G 转移到声子T_ph，当T_ph 达到熔点，例如2700 K，一些热量(360000 J) 潜伏融化前加热。

这是我的代码：

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 中添加一个高斯峰，并且峰中心应该在熔点附近。

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

请帮帮我

【问题讨论】：

请用数学说明这两种表示是什么
它们可能是在特定温度 T_ph 下的 delta 函数
请写下 math 明确你想要发生的事情。这不是解决问题物理问题的合适场所；这是一个研究问题。我将在这里提供帮助，将数学转换为 FiPy 语法。
抱歉耽搁了，我做了一些研究。我需要从 eq1 中减去这个函数：( L = (1/numerix.sqrt(2*numerix.pi * sig2)) * numerix.exp(-(T_ph - 1850)**2 / 2 * sig 2)) 总积分等于潜热(sig2 = 0.01)，它是一个高斯函数，但是它的峰很窄，怎么加？如何确保它被代码看到，使用的时间增量是多少？
请将其编辑到您的问题中。不保证评论会保留在 StackOverflow 上。

标签： python pde fipy

【解决方案1】：

根据 cmets（请将其编辑到问题中），将 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 更新时更新。

【讨论】：