Посылка
Я пытаюсь решить набор связанных 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
Как автор рассматриваемого программного обеспечения я не согласен. Этот вопрос представляет собой полностью проработанный пример и адекватное описание того, что они ищут. @ itprorh66: StackOverflow имеет заслуженную репутацию неприветливого человека; брось. — person nabz schedule 30.03.2021
Какую ошибку вы получаете? Что за бессмысленный результат? — person nabz schedule 31.03.2021
Это не было ошибкой. Решатель работал, но решение не имело физического смысла. — person nabz schedule 31.03.2021
Источник: