FiPy для потока заряженных частиц

Посылка

Я пытаюсь решить набор связанных PDE, который описывает диффузию заряженных частиц с различные коэффициенты диффузии с использованием FiPy. Конечная цель — получить профиль концентрации как для веществ, так и для электрического поля.

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

Заряженные частицы диффундируют от центра домена (левая граница) к стенкам домена (правая граница). Это соответствует граничному условию Дирихле (B.C.) на левой границе, где концентрация обоих видов = 0, и граничному условию Neumann B.C. на правой границе, где поток частиц равен 0 для описания радиальной симметрии. Поскольку заряженные частицы диффундируют с разной скоростью, возникает электрическое поле, возникающее из-за пространственного заряда. Электрическое поле ускоряет более медленные частицы и замедляет более быстрые частицы пропорционально величине поля.

P — это концентрация положительно заряженных частиц, а N — концентрация отрицательно заряженных частиц. E — электрическое поле пространственного заряда.

Проблема

Кажется, я не могу получить разумное решение из своего кода, и я думаю, что это может быть связано с тем, как я использовал условия градиента / расхождения как ConvectionTerm:

from fipy import *
import scipy.constants as constant
from fipy.tools import numerix
import numpy as np

## Defining physical constants
pi = constant.pi
m_argon = 6.6335e-26 # kg
k_b = constant.k # J/K
e_0 = constant.epsilon_0 # F/m
q_e = constant.elementary_charge # C
m_e = constant.electron_mass # kg
planck = constant.h

def char_diff_length(L,R):
    """Characteristic diffusion length in a cylinder.
    Used for determining the ambipolar diffusion coefficient.
    ref: https://doi.org/10.6028/jres.095.035"""
    a = (pi/L)**2
    b = (2.405/R)**2
    c = (a+b)**(1/2)
    return c

def L_Debye(ne,Te):
    """Electron Debye screening length given in m.
    ne is in #/m3, Te is in K."""
    if ne < 3.3e-5:
        ne = 3.3e-5
    return (((e_0*k_b*Te)/(ne*q_e**2)))**(1/2)

## Setting system parameters
# Operation parameters
Pressure = 1.e5 # ambient pressure Pa
T_g = 400. # background gas temperature K
n_g = Pressure/k_b/T_g # gas number density #/m3
Q_std = 300. # standard volumetric flowrate in sccm
T_e_0 = 11. # plasma temperature ratio T_e/T_g here assumed to be T_e = 0.5 eV and T_g = 500 K
n_e_0 = 1.e20 # electron density in bulk plasma #/m3

# Geometric parameters
R_b = 1.e-3 # radius cylinder m
L = 1.e-1 # length of cylinder m

# Transport parameters
D_ion = 4.16e-6 #m2/s ion diffusion, obtained from https://doi.org/10.1007/s12127-020-00258-z
mu_ion = D_ion*q_e/k_b/T_g # ion electrical mobility using Einstein relation

D_e = 100.68122*D_ion #m2/s electron diffusion
mu_e = D_e*q_e/k_b/T_g # electron electrical mobility using Einstein relation

Lambda = char_diff_length(L,R_b)
debyelength_e = L_Debye(n_e_0,T_g)
gamma = (Lambda/debyelength_e)**2
delta = D_ion/D_e



def d_j(rb,n): #sets the desired spatial steps for mesh
    dj = np.zeros(n)
    for j in range(n):
        dj[j] = 2*rb*(1 - j/n)/n
    return dj

#Initializing mesh
dj = d_j(1.,100) # 100 points
mesh = CylindricalGrid1D(dr = dj)

#Declaring cell variables
N = CellVariable(mesh=mesh, value = 1., hasOld = True, name = "electron density")
P = CellVariable(mesh=mesh, value = 1., hasOld = True, name = "ion density")
H = CellVariable(mesh=mesh, value = 0., hasOld = True, name = "electric field")

#Setting boundary conditions
N.constrain(0.,mesh.facesRight) # electron density = 0 at walls
P.constrain(0.,mesh.facesRight)# ion density = 0 at walls

H.constrain(0.,mesh.facesLeft) # electric field = 0 in the center
N.faceGrad.constrain([0.],mesh.facesLeft) # flux of electron = 0 in the center
P.faceGrad.constrain([0.],mesh.facesLeft) # flux of ion = 0 in the center

if __name__ == '__main__':
    viewer = Viewer(vars=(P,N))
    viewer.plot()

eqn1 = (TransientTerm(var=P) == DiffusionTerm(coeff=delta,var=P) 
                                - ConvectionTerm(coeff=[H.cellVolumeAverage,],var=P) 
                                - ConvectionTerm(coeff=[P.cellVolumeAverage,],var=H))
eqn2 = (TransientTerm(var=N) == DiffusionTerm(var=N) 
                                + (1/delta)*(ConvectionTerm(coeff=[H.cellVolumeAverage,],var=N)
                                            +ConvectionTerm(coeff=[N.cellVolumeAverage,],var=H)))
eqn3 = (TransientTerm(var=H) == gamma*(ConvectionTerm(coeff=[delta**2,],var=P) 
                                       - ConvectionTerm(coeff=[delta,],var=N) 
                                       - H*(delta*P.cellVolumeAverage + N.cellVolumeAverage)))
P.setValue(1.)
N.setValue(1.)
H.setValue(0.)
eqn1d = eqn1 & eqn2 & eqn3



timesteps = 1e-5
steps = 100
for i in range(steps):
    P.updateOld()
    N.updateOld()
    H.updateOld()
    res = 1e10
    sweep = 0
    while res > 1e-3 and sweep < 20:
        res = eqn1d.sweep(dt=timesteps)
        sweep += 1
    if __name__ == '__main__':
        viewer.plot()

добро пожаловать в Stack Overflow.! Ваш вопрос кажется слишком общим для этого сайта. Чтобы мы могли вам помочь, нам нужен воспроизводимый набор данных, минимальный код, который генерирует ошибку / проблему, что вы получаете и чего ожидаете. Прочтите С чего начать и Минимальный воспроизводимый пример, затем отредактируйте свое сообщение.   —  person nabz    schedule 30.03.2021

См. также:  Подпроцесс Python.Popen + ffmpeg прерывает ввод терминала

Как автор рассматриваемого программного обеспечения я не согласен. Этот вопрос представляет собой полностью проработанный пример и адекватное описание того, что они ищут. @ itprorh66: StackOverflow имеет заслуженную репутацию неприветливого человека; брось.   —  person nabz    schedule 30.03.2021

Какую ошибку вы получаете? Что за бессмысленный результат?   —  person nabz    schedule 31.03.2021

Это не было ошибкой. Решатель работал, но решение не имело физического смысла.   —  person nabz    schedule 31.03.2021

Понравилась статья? Поделиться с друзьями:
IT Шеф
Комментарии: 2
  1. nabz

    Электрическое поле — это вектор, а не скаляр.

    H = CellVariable(rank=1, mesh=mesh, value = 0., hasOld = True, name = "electric field")
    

    Исправление должно сделать преобразование терминов в FiPy более понятным:

    Нет причин запускать цепное правило на последнем члене eq1 или eq2; они уже находятся в канонической форме для FiPy ConvectionTerm. После цепного правила они становятся, например, , ни одна из этих форм не нравится FiPy. Вы можете написать эти два последних термина как явные источники, но не должны.

    eqn1 = (TransientTerm(var=P) == DiffusionTerm(coeff=delta, var=P) 
                                    - ConvectionTerm(coeff=H, var=P))
    eqn2 = (TransientTerm(var=N) == DiffusionTerm(var=N) 
                                    + (1/delta)*ConvectionTerm(coeff=H, var=N))
    

    Я не совсем понимаю eq3. Это похоже на интегрирование уравнения неразрывности? Я не вижу этого при быстром сканировании цитируемой вами статьи Фелпса. Как бы то ни было, FiPy не в такой форме; вы можете написать это, но это не решит хорошо. Термины справа не ConvectionTerms, это просто градиенты.

    Если вы собираетесь разрешить разделение зарядов и беспокоиться о длине Дебая, я думаю, вам следует решить уравнение Пуассона. Не могли бы вы рассказать, откуда взялось это уравнение? Возможно, мы сможем придать ему форму, которая понравится FiPy.

    Привет, Джонатан, спасибо за ответ. Я последую вашему совету, чтобы соблюдать правила до цепочки. person nabz; 31.03.2021

  2. nabz

    eq3 — модифицированное уравнение Пуассона. Я попытался выполнить процедуру, описанную Freeman, где производная по времени выполняется для уравнение Пуассона для замены уравнений неразрывности видов. Фриман решил эти уравнения, используя пакет Gear, который, как я могу предположить, является пакетом на Фортране. Я пошёл по его стопам из-за наивности, потому что не очень разбираюсь в численных методах.

    Я снова попробую решить с уравнением Пуассона в его стандартной форме.

    изменить: я изменил электрическое поле H на тензор ранга 1, и я изменил eq3, а также немного изменил определение gamma. Все остальное без изменений.

    H = CellVariable(rank = 1, mesh=mesh, value = 0., hasOld = True, name = "electric field")
    
    charlength = char_diff_length(L,R_b)
    debyelength_e = L_Debye(n_e_0,T_g)
    gamma = (debyelength_e/charlength)**2
    delta = D_ion/D_e
    
    eqn1 = (TransientTerm(var=P) == DiffusionTerm(coeff=delta,var=P) 
                                    - ConvectionTerm(coeff=H,var=P))
    eqn2 = (TransientTerm(var=N) == DiffusionTerm(var=N) 
                                    + (1/delta)*ConvectionTerm(coeff=H,var=N))
    eqn3 = (ConvectionTerm(coeff = gamma/delta, var=H) == ImplicitSourceTerm(var=P) 
                                    - ImplicitSourceTerm(var=N))
    P.setValue(1.)
    N.setValue(1.)
    H.setValue(0.)
    eqn1d = eqn1 & eqn2 & eqn3
    
    timesteps = 1e-8
    steps = 100
    for i in range(steps):
        P.updateOld()
        N.updateOld()
        H.updateOld()
        res = 1e10
        sweep = 0
        while res > 1e-3 and sweep < 20:
            res = eqn1d.sweep(dt=timesteps)
            sweep += 1
        if __name__ == '__main__':
            viewer.plot()
    

    Он не дает мне тех же ошибок, что и раньше, что является некоторым индикатором прогресса. Однако он выдает новую ошибку:

    ValueError: all the input arrays must have same number of dimensions, but the array at index 0 has 1 dimension(s) and the array at index 2 has 2 dimension(s)

    Я понимаю. Думаю, я вижу некоторую полезность в подходе Фримена, учитывая желание использовать решатель ODE почти пятьдесят лет назад. В наши дни решить уравнение Пуассона несложно. Более того, любая переходная версия уравнения Пуассона действует со скоростью света. Мне никогда не приходилось использовать что-либо, кроме квазистатической версии, кроме таких вещей, как полупроводниковые устройства ГГц или ТГц. person nabz; 31.03.2021

    Рекомендую решать уравнение Пуассона. Я рекомендую решать для потенциального (скалярного), а не электрического поля (вектора). $ \ vec {E} = — \ nabla V $. Это превращает уравнение Пуассона в форму УЧП 2-го порядка, что хорошо подходит для FiPy. person nabz; 31.03.2021

    Для дальнейшего обсуждения я рекомендую вам использовать наш список рассылки или средство отслеживания проблем GitHub. StackOverflow не очень подходит для постоянного обсуждения и устранения неполадок. person nabz; 31.03.2021

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

    спасибо, я буду следить за дальнейшим обсуждением в списке рассылки person nabz; 31.03.2021

Добавить комментарий

;-) :| :x :twisted: :smile: :shock: :sad: :roll: :razz: :oops: :o :mrgreen: :lol: :idea: :grin: :evil: :cry: :cool: :arrow: :???: :?: :!: