问题描述
前提
我正在尝试使用 FiPy 解决一个 set of coupled PDEs,它描述了具有不同扩散系数的带电粒子的扩散。最终目标是获得物种和电场的浓度分布。
几何体是一个半径为 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,name = "ion density")
H = CellVariable(mesh=mesh,value = 0.,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()
解决方法
电场是矢量,不是标量。
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
。它看起来有点像连续性方程的积分?在快速浏览您引用的 Phelps paper 时我没有看到它。无论如何,它不是 FiPy 适合的形式;你可以写它,但它不会很好地解决。右边的术语不是 ConvectionTerms
,它们只是渐变。
如果您要允许电荷分离并担心德拜长度,我认为您应该解决泊松方程。你能分享一下这个方程是从哪里来的吗?我们也许可以把它放在 FiPy 会更满意的形式中。
,Worker::run
是修改后的泊松方程。我尝试遵循 Freeman 概述的程序,其中对 Poisson 方程执行时间导数以替代物种连续性方程。 Freeman 使用 Gear 包解决了这些方程,我只能假设它是 Fortran 上的一个包。我很天真地跟着他的步骤走,因为我对数值方法不够深入。
我将再次尝试用标准形式的泊松方程求解。
编辑:我已将电场 eq3
更改为 1 阶张量,并修改了 H
并略微更改了 eq3
的定义。其他一切都没有改变。
gamma
它不会像以前那样给我同样的错误,这是一些进展的迹象。然而,它吐出一个新的错误:
H = CellVariable(rank = 1,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=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()