キャリア連続の式をpythonでシミュレーションしたい

Question

キャリア連続の式をpythonでシミュレーションしたい

pythonを用いてキャリア連続の式（電流連続の式）をシミュレーションしています。
ポアソン方程式と電界の式とキャリア連続の式はそれぞれに相互関係があり、
その前段階としてポアソン方程式と電界の式をシミュレーションをしました。
ポアソン方程式と電界の式のシミュレーションはネットにあるものを参考にして成功したのですが、
キャリア連続の式のシミュレーション（離散化や数値の設定）について参考文献がなく
シミュレーションができませんでした。

発生している問題・エラーメッセージ

一面同一色の画像が得られる。

該当のソースコード

k_b = 1.38*10**-23 # ボルツマン定数
e = 1.60*10**-19 # 電気素量
T = 300 # 絶対温度 
myu_e = 100 # 移動度
n_0 = 10**-3 # 平衡時の正孔の密度 

# 拡散係数D_eを求める
D_e = float(myu_e*(k_b*T / e))

tau = 10**-4 # 少数キャリアの寿命
convegence_criterion = 10**-5

# キャリア連続の式の計算
n = np.zeros([Lx + 1, Ly + 1])
ni_x = np.zeros([Lx + 1, Ly + 1])
ni_y = np.zeros([Lx + 1, Ly + 1])

def initialize_n(condition, n, delta_Lx, delta_Ly, Lx, Ly):
    for i in range(Lx):
        for j in range(Ly):
            if condition(i, j, delta_Lx, delta_Ly):
                n[i, j] = 1

def evolve_ts(im, dx2, dy2, Lx, Ly, myu_e, Ex, Ey, D_e, tau, n_0, dt, delta_Lx):
    global n, ni_x
    n[1:-1, 1:-1] = ni_x[1:-1, 1:-1] + dt*(myu_e*Ex*((ni_x[2:, 1:-1] - ni_x[-2:, 1:-1])/(2*delta_Lx)) + 
                                           D_e*((ni_x[2:, 1:-1] - 2*ni_x[1:-1, 1:-1] + ni_x[-2:, 1:-1])/dx2) - 
                                         (ni_x[1:-1, 1:-1] - n_0)/tau)
    ni_x = np.copy(n)
    im.set_array(ni_x)
    return im, 

def main():
    global n, ni_x, delta_Lx, delta_Ly, Lx, Ly
    dx2 = delta_Lx ** 2 
    dy2 = delta_Ly ** 2

    dt = dx2 * dy2 / (2 * (dx2 + dy2))  

    condition = (lambda i, j, delta_Lx, delta_Ly:
                 (((i*delta_Lx - 0.5)**2 + (j*delta_Ly - 0.5)**2 <= 0.01)
                 & ((i*delta_Lx - 0.5)**2 + (j*delta_Ly - 0.5)**2 >= .1))
                 )

    initialize_n(condition, ni_x, delta_Lx, delta_Ly, Lx, Ly)

    fig = plt.figure()
    img = plt.subplot(111)
    im = img.imshow(ni_x, cmap = 'hsv', interpolation = 'nearest', origin = 'lower')
    im.figure = fig
    fig.colorbar(im)
    ani = animation.FuncAnimation(fig, evolve_ts, 240, 
                                  fargs = (im, dx2, dy2, Lx, Ly, dt, D_e, tau,),
                                  blit = True, interval = 50)
    
    plt.show()

if __name__ == "__main__":
    main()



# Commented out IPython magic to ensure Python compatibility.
# %%time
#

試したこと

conditionや移動度、平衡時の正孔の密度、少数キャリアの寿命の値を変えた。

補足情報（FW/ツールのバージョンなど）

ここにより詳細な情報を記載してください。

Comments

[jbpb0]

物理現象の理論式が正しくコードに実装されてるかみたいな話だと思うので、理論式が分かる情報が書かれてるWebページのURLとかもあると、回答されやすくなるように思います

[nagi-kimura]

(dn/dt)=μEx(dn/dx)+De(d^2n/dx^2)-(n-n0)/τe
(dn/dt)=μEy(dn/dy)+De(d^2n/dy^2)-(n-n0)/τe
dは全てラウンドd、μは移動度、Exは電場のx成分、Eyは電場のy成分、nは正孔密度、Deは拡散係数、n0は平衡時における正孔密度、τeは正孔の寿命です。
正孔密度の時間tにおける二階偏微分方程式のため見づらくなり申し訳ありません。

[jbpb0]

> キャリア連続の式のシミュレーション（離散化や数値の設定）について参考文献がなく

https://core.ac.uk/download/pdf/267937503.pdf
の式(26)が近いのかな?
違うかな?

[nagi-kimura]

拝読しました。
私が得たい離散化した数式に酷似してました。
しかし少し数式が異なることから、多分正孔と電子のキャリア密度を同時に表した物ではないかと思いました。
大変貴重な文献をいただき心から感謝申し上げます。

Answer 1

evolve_ts() の2行目の計算は、理論式のx方向の式
(dn/dt)=μEx(dn/dx)+De(d^2n/dx^2)-(n-n0)/τe
に相当するのだと思うのですが、比べてて変だなと思うところがあります

たとえば、μEx(dn/dx) に相当する

myu_e*Ex*((ni_x[2:, 1:-1] - ni_x[-2:, 1:-1])/(2*delta_Lx))

ですが、ni_x の最初の次元の添え字が 2: だと、最初の二つの要素以外の全部、という意味になります
一方、-2: だと、最後の二つの要素だけ、という意味になります

pythonで下記を実行して、確認してみてください

import numpy as np
xxx = np.arange(10)
print(xxx)
print(xxx[2:])
print(xxx[-2:])

微分を中心差分で代用しているのが意図だと思うので、二つ目の ni_x の添え字は -2: ではなく :-2 ではないですか?

print(xxx[2:])
print(xxx[:-2])

De(d^2n/dx^2) のところも、三つ目の ni_x の添え字が同様になってます

Comments

[nagi-kimura]

変だなと思わせて申し訳ありません。まだPythonを始めて1ヶ月程しか経ってなく離散化についても講義でほぼ触れられなかったため、離散化についてもインターネットの情報をもとにして行いました。訂正していただきありがとうございます。