FiPy を使用して太陽電池をシミュレートしようとしていますが、単純なテスト ケースでも妥当な結果を得るのに苦労しています。
私のテストの問題は、平衡状態の暗闇での突然の 1D pn ホモ接合です。方程式の支配系は、追加の生成または再結合のない半導体方程式です。
ポアソン方程式は、電子 ( n )、正孔 ( p )、ドナー ( N D )、およびアクセプター ( N A )の密度が与えられると、誘電率εを持つ半導体の電界 ( φ ) を決定します。電子はq :
∇² φ = q ( p − n + N D − N A ) / ε
電子と正孔は、移動度 ( μ ) と拡散定数 ( D )に応じて、電流密度Jでドリフトおよび拡散します。
J n = qμ n n E + qD n ∇n
J p = qμ p p E − 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)
メッシュ サイズ、時間ステップ、時間ステップ数、スイープ数などのさまざまな値を試してみました。多少の変動は見られますが、現実的な解決策を提供する一連の条件を見つけることができませんでした。問題はおそらく現在の用語の表現にあると思います。
通常、これらの方程式を解く場合、電流密度は、直接離散化ではなく、Scharfetter-Gummel (SG) 離散化スキームを使用して近似されます。SG 方式では、セル面を通過する電子電流密度 ( J ) は、セルKおよびLのいずれかの側の中心で定義された電位 ( φ ) および電荷密度 ( n )の値の関数として近似されます。
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 を使用して半導体デバイスをシミュレートできる代替アプローチはありますか?