时间戳

首页 > 解决方案 > 如何在 FiPy 中添加潜热项？

python - 如何在 FiPy 中添加潜热项？

问题描述

我一直在尝试使用模型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传递到声子的热量，例如，当达到熔点时，一些热量在熔化前作为潜热。T_phGT_ph2700 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)/√22
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峰中心应该在熔点附近。

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

请帮帮我

标签： pythonpdefipy

解决方案

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

推荐阅读