Как решить уравнения, используя схему Шарфеттера-Гуммеля в FiPy?

Question

Я пытаюсь использовать FiPy для моделирования солнечных элементов, но я изо всех сил пытаюсь получить разумные результаты даже для простых тестовых случаев.

Моя тестовая задача - резкий 1D pn гомопереход в темноте в равновесии. Управляющей системой уравнений являются уравнения полупроводников без дополнительной генерации или рекомбинации.

Уравнение Пуассона определяет электрическое поле (φ) в полупроводнике с диэлектрической проницаемостью & epsi;, учитывая плотности электронов (n), дырок (p), доноров (N_D) и акцепторов (N_A), где заряд электрона равен q:

∇²φ = q(p - n + N_D - N_A) / & epsi;

Электроны и дырки смещаются и диффундируют с плотностями тока, J, в зависимости от их подвижностей (μ) и констант диффузии (D):

J_n = qμ_n nE + qD_n∇n

J_p = qμ_p pE - qD_p∇n

Эволюция заряда в системе учитывается уравнениями неразрывности электронов и дырок:

∂n / ∂t = (∇ ·J_n) / q

∂p / ∂t = - (∇ ·J_p) / q

которая может быть выражена в канонической форме FiPy как:

∂n / ∂t = μ_n∇ · (-n∇φ) + D_n∇²n

∂p / ∂t = - (μ_p∇ · (-p∇φ) - D_p∇²n)

Чтобы попытаться решить проблему в FiPy, я сначала импортирую модули и определяю физические параметры.

from __future__ import print_function, division
import fipy
import numpy as np
import matplotlib.pyplot as plt

eps_0 = 8.8542e-12     # Permittivity of free space, F/m
q     = 1.6022e-19     # Charge of an electron, C
k     = 1.3807e-23     # Boltzmann constant, J/K
T     = 300            # Temperature, K
Vth   = (k*T)/q        # Thermal voltage, eV

N_ap  = 1e22           # Acceptor density in p-type layer, m^-3
N_an  = 0              # Acceptor density in n-type layer, m^-3
N_dp  = 0              # Donor density in p-type layer, m^-3
N_dn  = 1e22           # Donor density in n-type layer, m^-3

mu_n  = 1400.0e-4      # Mobilty of electrons, m^2/Vs
mu_p  = 450.0e-4       # Mobilty of holes, m^2/Vs
D_p   = k*T*mu_p/q     # Hole diffusion constant, m^2/s
D_n   = k*T*mu_n/q     # Electron diffusion constant, m^2/s
eps_r = 11.8           # Relative dielectric constant

n_i   = (5.29e19 * (T/300)**2.54 * np.exp(-6726/T))*1e6       
V_bias = 0             

Затем создайте сетку, переменные решения и профиль легирования.

nx = 20000
dx = 0.1e-9
mesh = fipy.Grid1D(dx=dx, nx=nx)
Ln = Lp = (nx/2) * dx

phi = fipy.CellVariable(mesh=mesh, hasOld=True, name='phi')
n = fipy.CellVariable(mesh=mesh, hasOld=True, name='n')
p = fipy.CellVariable(mesh=mesh, hasOld=True, name='p')
Na = fipy.CellVariable(mesh=mesh, name='Na')
Nd = fipy.CellVariable(mesh=mesh, name='Nd')

Затем я устанавливаю некоторые начальные значения для центров клеток и накладываю граничные условия Дирихле на все параметры.

x = mesh.cellCenters[0]

n0 = n_i**2 / N_ap
nL = N_dn
p0 = N_ap
pL = n_i**2 / N_dn
phi_min = -(Vth)*np.log(p0/n_i)
phi_max = (Vth)*np.log(nL/n_i) + V_bias

Na.setValue(N_an, where=(x >= Lp))
Na.setValue(N_ap, where=(x < Lp))
Nd.setValue(N_dn, where=(x >= Lp))
Nd.setValue(N_dp, where=(x < Lp))
n.setValue(N_dn, where=(x > Lp))
n.setValue(n_i**2 / N_ap, where=(x < Lp))
p.setValue(n_i**2 / N_dn, where=(x >= Lp))
p.setValue(N_ap, where=(x < Lp))
phi.setValue((phi_max - phi_min)*x/((Ln + Lp)) + phi_min)

phi.constrain(phi_min, mesh.facesLeft)
phi.constrain(phi_max, mesh.facesRight)
n.constrain(nL, mesh.facesRight)
n.constrain(n_i**2 / p0, mesh.facesLeft)
p.constrain(n_i**2 / nL, mesh.facesRight)
p.constrain(p0, mesh.facesLeft)

Я выражаю уравнение Пуассона как

eps = eps_0*eps_r
rho = q * (p - n + Nd - Na)
rho.name = 'rho'

poisson = fipy.ImplicitDiffusionTerm(coeff=eps, var=phi) ==  -rho

уравнения неразрывности как

cont_eqn_n = (fipy.TransientTerm(var=n) == 
              (fipy.ExponentialConvectionTerm(coeff=-phi.faceGrad*mu_n, var=n) 
               + fipy.ImplicitDiffusionTerm(coeff=D_n, var=n)))

cont_eqn_p = (fipy.TransientTerm(var=p) == 
              - (fipy.ExponentialConvectionTerm(coeff=-phi.faceGrad*mu_p, var=p) 
               - fipy.ImplicitDiffusionTerm(coeff=D_p, var=p)))

и решить, связав уравнения и развертки:

eqn  = poisson & cont_eqn_n & cont_eqn_p

dt = 1e-12
steps = 50
sweeps = 10    

for step in range(steps):
    phi.updateOld()
    n.updateOld()
    p.updateOld()
    for sweep in range(sweeps):
        eqn.sweep(dt=dt)

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

Обычно при решении этих уравнений плотности тока аппроксимируются с использованием схемы дискретизации Шарфеттера-Гуммеля (SG), а не прямой дискретизации. В схеме SG плотность электронного тока (J) через поверхность ячейки аппроксимируется как функция значений потенциала (φ) и плотности заряда (n), определенных на центрах ячеек K и L с обеих сторон, как

J_{n, KL}=qμ_n V_T[B(δφ / V_T)n_L - B(-δφ / V_T)n_K)

где q - заряд электрона, µ_n - подвижность электрона, V_T - тепловое напряжение, δφ=φ_L-φ_K, а B(x) - функция Бернулли x/ (e^x−1).

Для меня не очевидно, как реализовать схему в FiPy. Я видел, что есть scharfetterGummelFaceVariable но я не могу понять из документации, подходит ли она или предназначена для этой проблемы. Глядя на код, кажется, что он вычисляет только функцию Бернулли, умноженную на коэффициент e^φ_L. Можно ли напрямую использовать scharfetterGummelFaceVariable решить этот тип проблемы? Если так, то как? Если нет, есть ли альтернативный подход, который позволит мне моделировать полупроводниковые устройства с использованием FiPy?