Блог

Главная - Вопросы по программированию - Как добавить скрытый термин тепла в FiPy?

Как добавить скрытый термин тепла в FiPy?

#python #pde #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_ph G , 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 amp; 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)

Как я понял, я могу включить скрытое тепло в качестве термина источника в качестве стока в уравнении 1 или добавить гауссов пик C_ph , и центр пика должен быть около точки плавления.

Я понятия не имею, какой из них лучше и стабильнее, я понятия не имею, как реализовать любой из них .

пожалуйста, помогите мне с этим

Комментарии:

1. Пожалуйста, покажите с помощью математики, что представляют собой эти два представления

2. они могут быть дельта-функцией, действующей при определенной температуре T_ph

3. Пожалуйста, запишите математику для того, что вы хотите, чтобы произошло. Это неподходящее место для того, чтобы разгадывать физику вашей проблемы; это исследовательская проблема. Я предложу здесь помощь в переводе математики на синтаксис FiPy.

4. Извините за задержку, я провел кое-какие исследования. Мне нужно вычесть эту функцию из eq1:( Л = (1/numerix.функция sqrt(2*numerix.Пи * сиг 2)) * numerix.ехр(-(T_ph — 1850)**2 / 2 * сиг 2)) с полной интеграции равна скрытой теплоте(sig2 = 0.01), его функцией Гаусса, но ее узкий пик, как добавить его? как убедиться, что это было замечено кодом с любым используемым приращением времени?

5. Пожалуйста, отредактируйте это в своем вопросе. Комментарии не гарантированно будут сохранены в StackOverflow.

Ответ №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 каждом обновлении.

Метки: