Computational algorithms for a one-dimensional MEP model of charge transfer in semiconductors

Full Text

Введение

В настоящий момент существует огромное число математических моделей, предназначенных для описания физических явлений в полупроводниковых устройствах. Поскольку поиск приближенных решений затруднителен, то встает вопрос о разработке численных алгоритмов нахождения приближенных решений. В работе рассматривается одномерный вариант MEP (Maximum Entropy Principle) модели [1].

Безразмерная математическая модель, описывающая перенос заряда в полупроводнике имеет следующий вид [2]:

$\left.\begin{array}{l}{R}_t+J_x\text{ =0,}\\ J_t+{\left(\frac{\text{2}}{\text{3}}RE\right)}_x\text{=}RQ\text{+}{c}_{11}J+c_{12}I,\\ {\left(RE\right)}_t+I_x\text{ =}JQ\text{+}cP\text{,}\\ I_t+{\left(\frac{\text{10}}{\text{9}}R{E}^2\right)}_x\text{=}\frac{\text{5}}{\text{3}}REQ\text{+}{c}_{21}J+c_{22}I,\end{array}\right\}$                           (1)

$\epsilon {\phi }_{xx}=R-\rho .$                                                                                          (2)

Где R – плотность электронов, u – скорость электронов, q – поток энергии, J = Ru, I = Rq, φ – электрический потенциал, Q = φ_x, E – энергия электронов, $\sigma =\frac{2}{3}E-1$, $P=R\sigma$, ε > 0 – константа, ρ = ρ(x), c, c₁₁, c₁₂, c₂₁, c₂₂ – коэффициенты, зависящие от энергии E [2, 3].

Системе уравнений (1) задали граничные условия (3), начальные данные (4), а уравнению Пуассона (2) поставили краевые условия (5):

$\left.\begin{array}{l}R\left(t\mathrm{,0}\right)=R\left(t\mathrm{,1}\right)=\mathrm{1,}\\ E\left(t\mathrm{,0}\right)=E\left(t\mathrm{,1}\right)=\frac{3}{2},\end{array}\right\}$                                                                             (3)

$\left.\begin{array}{l}R\left(\mathrm{0,}x\right)=R{}_0\left(_x^{\right)},\\ J\left(\mathrm{0,}x\right)=J{}_0\left(_x^{\right)},\\ E\left(\mathrm{0,}x\right)=E{}_0\left(_x^{\right)},\\ I\left(\mathrm{0,}x\right)=I{}_0\left(_x^{\right)},\end{array}\right\}$                                                                                    (4)

$\phi \left(t\mathrm{,0}\right)=\mathrm{0,}\phi \left(t\mathrm{,1}\right)=B,$                                                                              (5)

$B>\mathrm{0,} R{}_0\left(x\right)>\mathrm{0,}E{}_0\left(x\right)>0$    .                       

Сформулированная задача (1) – (5) является одномерной задачей о баллистическом диоде. Это полупроводниковый прибор, состоящий из трех областей. Первую и третью области называют n⁺ - зонами. Они имеют высокую концентрацию легирования. Средняя область имеет низкую концентрацию легирования и называется n - зоной.

Безразмерная величина ρ имеет следующий вид:

$\rho =\left\{\begin{array}{l}\mathrm{1,}\text{если }x\text{ принадлежи }{n}^{+}-\text{зоне,}\\ \delta =\frac{N^{+}}{N},\text{если }x\text{ принадлежи }n-\text{зоне,}\end{array}\right.$                                                                             

N, N⁺ являются плотностями легирования в n - зоне, n⁺ -зоне. График ρ(x) изображен на рисунке 1.

Рисунок 1 – Схематическое изображение n⁺– n – n⁺ баллистического диода

При осуществлении поиска приближенных решений в данной задаче возникли вычислительные проблемы, связанные с тем, что уравнения модели являются нелинейными уравнениями и содержат очень малые параметры. В задаче не удалось найти приближенные решения с помощью разностных схем. Основные проблемы возникли с нахождением энергии электронов. Энергия электронов очень сильно увеличивалась, и происходило переполнение системы.

Результаты и обсуждение

Для нахождения приближенных решений задачу (1) – (5) в стационарном случае удалось свести к системе из трех уравнений Пуассона [4]:

$\left\{\begin{array}{l}\frac{d^2\phi }{d{x}^2}={\mathcal{F}}^{\left(\phi \right)}(\chi ,\rho )=\beta (e^{\chi }-\rho ),\\ \frac{d^2\sigma }{d{x}^2}={\mathcal{F}}^{\left(\sigma \right)}\left(\Sigma ,X,Q,\sigma \right)=a_1{\Sigma }^2+a_2\Sigma X+a_3\Sigma Q+a_{4}XQ+a_5{Q}^2+bc\sigma ,\\ \frac{d^2\chi }{d{x}^2}={\mathcal{F}}^{\left(\chi \right)}\left(\Sigma ,X,Q,\sigma ,X,\rho \right)=-X^2+b_1{\Sigma }^2+b_2\Sigma X+b_3\Sigma Q+b_{4}XQ+\\ +b_5{Q}^2+\frac{\beta }{1+\sigma }\left(e^{\chi }-\rho \right)+nc\sigma ,\end{array}\right.$         (6)

$\left\{\begin{array}{l}\chi =\sigma =0\text{ при }x=\mathrm{0,1,}\\ \phi =0\text{ при }x=\mathrm{0,}\\ \phi =B\text{ при }x=\mathrm{1,}\end{array}\right.$                                                                         (7)

$\begin{array}{l}u=F\left(E\right)\left\{Q-\left(1+\sigma \right)X-F_0\left(E\right)\Sigma \right\},\\ q=G\left(E\right)\left\{-Q+\left(1+\sigma \right)X+G_0\left(E\right)\Sigma \right\}.\end{array}$                                                         (8)

Здесь $\beta =\frac{1}{\epsilon }, \chi =\ln R, \Sigma =\frac{d\sigma }{dx}, X=\frac{d\chi }{dx}$

$\begin{array}{c}{a}_1=-a^{\text{'}}F\left(E\right)F_0\left(E\right)+b^{\text{'}}G\left(E\right)G_0\left(E\right),\\ a_2=-1+\left(1+\sigma \right)\left\{b^{\text{'}}G\left(E\right)-a^{\text{'}}F\left(E\right)\right\},\\ a_3=a^{\text{'}}F\left(E\right)-b^{\text{'}}G\left(E\right)-bF\left(E\right)F_0\left(E\right),\\ a_4=-b\left(1+\sigma \right)F\left(E\right),\\ a_5=bF\left(E\right)\end{array}$

$\begin{array}{c}{b}_1=-m^{\text{'}}F\left(E\right)F_0\left(E\right)+n^{\text{'}}G\left(E\right)G_0\left(E\right),\\ b_2=\left(1+\sigma \right)\left\{n^{\text{'}}G\left(E\right)-m^{\text{'}}F\left(E\right)\right\},\\ b_3=-\frac{1}{{\left(1+\sigma \right)}^2}+m^{\text{'}}F\left(E\right)-n^{\text{'}}G\left(E\right)-nF\left(E\right)F_0\left(E\right),\\ b_4=\frac{1}{1+\sigma }+bF\left(E\right)F_0\left(E\right),\\ b_5=nF\left(E\right)\end{array}$

$\begin{array}{c}a=a\left(E\right)=\frac{2}{5}\frac{c_{21}}{1+\sigma }-c_{11},\\ b=b\left(E\right)=\frac{2}{5}\frac{c_{22}}{1+\sigma }-c_{12},\\ m=m\left(E\right)=\frac{c_{11}-a}{1+\sigma },\\ n=n\left(E\right)=\frac{c_{12}-b}{1+\sigma },\end{array}$

$\begin{array}{c}F\left(E\right)=-\frac{c_{22}-\frac{5}{3}E{c}_{12}}{det}, G\left(E\right)=-\frac{c_{21}-\frac{5}{3}E{c}_{11}}{det},\\ F_0\left(E\right)=1-\frac{\frac{5}{3}E{c}_{12}}{c_{22}-\frac{5}{3}E{c}_{12}}, G_0\left(E\right)=1-\frac{\frac{5}{3}E{c}_{11}}{c_{21}-\frac{5}{3}E{c}_{11}}, \\ det=c_{11}{c}_{22}-c_{21}{c}_{12} , a^{\text{'}}=\frac{da}{d\sigma }=\frac{3}{2}\frac{da}{dE} \mathrm{и} \mathrm{т}.\mathrm{д}.\end{array}$

Для поиска стационарных решений задачи (1) – (5) необходимо осуществить:

1. Задание начальных данных:

$\phi , E, \sigma =\frac{2}{3}E-1$.

2. Расчет необходимых коэффициентов:

$\begin{array}{c}c,c_{11},c_{12},c_{21},c_{22}, K\left(x\right)=\underset{0}{\overset{x}{\int }}\rho c\sigma d\tau , a_{11}\left(x\right)=\underset{0}{\overset{x}{\int }}\frac{c_{11}+c_{12}\phi }{\rho }d\tau ,\\ a_{12}\left(x\right)=\underset{0}{\overset{x}{\int }}\frac{c_{12}}{\rho }d\tau , a_{21}\left(x\right)=\underset{0}{\overset{x}{\int }}\frac{c_{21}+c_{22}\phi }{\rho \left(1+\sigma \right)}d\tau ,\\ a_{22}\left(x\right)=\underset{0}{\overset{x}{\int }}\frac{c_{22}}{\rho \left(1+\sigma \right)}d\tau , d_1\left(x\right)=\underset{0}{\overset{x}{\int }}\left\{\left(1+\sigma \right)\frac{{\rho }^{\text{'}}}{\rho }-\frac{c_{12}}{\rho }K\right\}d\tau ,\\ d_2\left(x\right)=\underset{0}{\overset{x}{\int }}\left\{\frac{5}{2}\left(1+\sigma \right)\frac{{\rho }^{\text{'}}}{\rho }-\frac{c_{22}}{\rho \left(1+\sigma \right)}K\right\}d\tau ,\\ \det =a_{11}\left(1\right)a_{22}\left(1\right)-a_{12}\left(1\right)a_{21}\left(1\right).\end{array}$

3. Вычисление промежуточных величин:

$\begin{array}{c}J=\frac{1}{\det }\left\{a_{22}\left(1\right)\left(d_1\left(1\right)-B\right)-a_{12}\left(1\right)\left(d_2\left(1\right)-\frac{5}{2}B\right)\right\},\\ I_0=\frac{1}{\det }\left\{-a_{21}\left(1\right)\left(d_1\left(1\right)-B\right)+a_{11}\left(1\right)\left(d_2\left(1\right)-\frac{5}{2}B\right)\right\}.\end{array}$

4. Пересчет φ , σ:

$\begin{array}{c}\omega \left(x\right)=d_1\left(x\right)-a_{11}\left(x\right)J-a_{12}\left(x\right)I_0^{},\\ \sigma \left(x\right)=\omega \left(x\right)-\frac{2}{5}{d}_2\left(x\right)+\frac{2}{5}{a}_{21}\left(x\right)J+\frac{2}{5}{a}_{22}\left(x\right)I_0^{},\\ \phi \left(x\right)=\omega \left(x\right)+\sigma \left(x\right), E\left(x\right)=\frac{3\left(\sigma \left(x\right)+1\right)}{2}.\end{array}$

В цикле организуется вычисление шагов 2−4. Цикл прекращает работу в тот момент, когда значения неизвестных найдены с заданной точностью.

Поиск приближенного решения задачи (6)−(7) осуществлялся с помощью метода установления [5]. Для этого рассматривали параболическую, гиперболическую и Соболевскую регуляризации [6, 7]. Было проведено огромное число экспериментов, которые показали, что скорость сходимости метода зависит от выбора регуляризации. Поскольку Соболевская регуляризация дает самую высокую скорость сходимости, то для построения численных алгоритмов используем следующую систему:

$\left\{\begin{array}{l}\left(1-{\xi }^2\right)\tau \phi ={\xi }^2\phi -{\mathcal{F}}^{\left(\phi \right)}(\chi ,\rho ),\\ \left(1-{\xi }^2\right)\tau \sigma ={\xi }^2\sigma -{\mathcal{F}}^{\left(\sigma \right)}(\xi \sigma ,\xi \chi ,\xi \phi ,\sigma ),\\ \left(1-{\xi }^2\right)\tau \chi ={\xi }^2\chi -{\mathcal{F}}^{\left(\chi \right)}(\xi \sigma ,\xi \chi ,\xi \phi ,\sigma ,\chi ,\rho ),\end{array}\right.$                                               (9)

$\tau =\frac{\partial }{\partial t}, \xi =\frac{\partial }{\partial x}, t -$ переменная по времени.

Начальные данные для систем уравнений (7), (9) при t > 0 и 0 < x < 1 примут следующий вид:

$\left\{\begin{array}{l}\phi \left(\mathrm{0,}x\right)={\phi }_0\left(x\right),\\ \sigma \left(\mathrm{0,}x\right)={\sigma }_0\left(x\right)=\frac{2}{3}{E}_0\left(x\right)-\mathrm{1,}\\ \chi \left(\mathrm{0,}x\right)={\chi }_0\left(x\right)=\ln R_0\left(x\right).\end{array}\right.$                                                           (10)

Далее для поиска численного решения задачи (7), (9) − (10) были сконструированы алгоритмы. В этих алгоритмах использовались кубические сплайны класса C² [8], интегральные уравнения и схема предиктор – корректор. Разработанные алгоритмы будем описывать на модельной задаче.

Система уравнений (9) состоит из трех одинаковых уравнений, которые отличаются друг от друга только правыми частями. Поэтому далее рассматриваем только одно из трех уравнений:

$\left\{\begin{array}{l}\left(1-{\xi }^2\right)\tau y={\xi }^{2}y-f^{\left(y\right)},t>\mathrm{0,}0<x<1;\\ y=0\text{ при }x=\mathrm{0,1,}t>0;\\ y=y_0\left(x\right)\text{ при }t=\mathrm{0,}0<x<1.\end{array}\right.$                                                 (11)

При подстановке функций ${\mathcal{F}}^{\left(\phi \right)}$, ${\mathcal{F}}^{\left(\sigma \right)}$, ${\mathcal{F}}^{\left(\chi \right)}$ в правую часть $f^{\left(y\right)}$, а при подстановке φ, σ, χ в функцию y, будет получена исходная система (9).

Также нам понадобиться и другая запись полученной задачи (11), которая нужна для сведения этой задачи к интегральным уравнениям:

$\left\{\begin{array}{l}\tau {\theta }_y+{\theta }_y={\tilde{f}}^{\left(y\right)},t>\mathrm{0,}0<x<1;\\ {\theta }_y={\theta }_{y_0}\left(x\right)\text{=}{y}_0\left(x\right)-{\xi }^2{y}_0\left(x\right)\text{ при }t=\mathrm{0,}0<x<\mathrm{1,}\end{array}\right.$                                   (12)

$\begin{array}{c}{\theta }_y={\theta }_y\left(t,x\right)\text{=}y\left(t,x\right)-{\xi }^{2}y\left(t,x\right),\\ {\tilde{f}}^{\left(y\right)}={\tilde{f}}^{\left(y\right)}\left(t,x\right)\text{=}y\left(t,x\right)-{\tilde{f}}^{\left(y\right)}\left(t,x\right).\end{array}$

Далее в уравнении $\left(1-{\xi }^2\right)\tau y={\xi }^{2}y-f^{\left(y\right)}$ системы (11) провели дискретизацию по переменной t. Производную τ y аппроксимировали выражением $\frac{Y-y}{\Delta }$:

${\xi }^{2}Y=\mathfrak{B}Y+\mathfrak{F},\mathfrak{B}=\frac{1}{1+\Delta },\mathfrak{F}=\frac{{\xi }^{2}y-y+\Delta f^{\left(y\right)}}{1+\Delta },$                                        (13)

 $∆$является шагом разностной сетки по времени t,

$y^n\left(x\right)=y\left(n\Delta ,x\right)=y,n=\mathrm{0,1,...,}Y=y^{n+1}\left(x\right)$.

Приближённое решение задачи (11) искали в виде кубического интерполяционного сплайна класса $C^2$(см. [8]):

$S\left(x\right)=(1-\overline{\tau })Y_k+\overline{\tau }{Y}_{k+1}-\frac{h_{}^2}{6}\overline{\tau }(1-\overline{\tau })\left[(2-\overline{\tau })m_k+(1+\overline{\tau })m_{k+1}\right],$                  (14)

$\begin{array}{c}\overline{\tau }=\frac{x-x_k}{h},x\in \left[x_k,x_{k+1}\right],k=\overline{\mathrm{0,}N-1}, \\ x_k = k . h, N . h = 1, Y_k = Y \left(x_k\right), m_k={\xi }^{2}Y\left(x_k\right).\end{array}$

Поскольку

$\xi S\left(x\right)=\frac{Y_{k+1}-Y_k}{h}-\frac{h_{}^{}}{6}\left[(2-6\overline{\tau }+3{\overline{\tau }}^2)m_k+(1-3{\overline{\tau }}^2)m_{k+1}\right],$

$x\in \left[x_k,x_{k+1}\right]$, $k=\overline{\mathrm{0,}N-1}$, то вычисляя

$\begin{array}{c}\xi S(x_k+0)=\frac{Y_{k+1}-Y_k}{h}-\frac{h_{}^{}}{6}\left[2m_k+m_{k+1}\right],\\ \xi S(x_k-0)=\frac{Y_k-Y_{k-1}}{h}-\frac{h_{}^{}}{6}\left[m_{k-1}+2m_k\right]\end{array}$

и приравнивая их, получили

$\frac{1}{2}{m}_{k-1}+2m_k+\frac{1}{2}{m}_{k+1}=\frac{3}{h^2}\left(\eta Y_k-\overline{\eta }{Y}_k\right), k=\overline{\mathrm{1,}N-1},$  ,                         (15)

$\eta =\psi -1, \overline{\eta }=1-{\psi }^{-1}, {\psi }^{\pm 1}{Y}_k=Y_{k\pm 1}\left({\psi }^{+1}=\psi \right)$

являются разностными операторами и операторами сдвига, соответственно.

Полагая в (13) x = x_k и подставляя ${\xi }^2{Y}_k$ в (15), переходим к следующей разностной схеме, которая является трехслойной:

$\begin{array}{l}\left\{1-\frac{h_{}^2}{6}\mathfrak{B}\right\}{Y}_{k-1}-2\left\{1+\frac{h_{}^2}{3}\mathfrak{B}\right\}{Y}_k+\left\{1-\frac{h_{}^2}{6}\mathfrak{B}\right\}{Y}_{k+1}=\\ =\frac{h_{}^2}{6}\left\{{\mathfrak{F}}_{k-1}^{}+4{\mathfrak{F}}_k+{\mathfrak{F}}_{k+1}\right\},k=\overline{\mathrm{1,}N-1},{\mathfrak{F}}_k=\mathfrak{F}\left(x_k\right).\end{array}$                             (16)

Для (16) задаем краевые условия:

$Y_0=\mathrm{0,}{Y}_N=0$.                                                                                      (17)

Решая трехдиагональную систему уравнений (16) − (17) с помощью метода прогонки [5], находим функции $Y_k,k=\overline{\mathrm{0,}N}$.

Для вычисления производной $\xi y$ также использовались кубические сплайны (14). Из

$\xi S\left(x\right)=\frac{Y_{k+1}-Y_k}{h}-\frac{h_{}^{}}{6}\left[(2-6\overline{\tau }+3{\overline{\tau }}^2)m_k+(1-3{\overline{\tau }}^2)m_{k+1}\right],$

$x\in \left[x_k,x_{k+1}\right], k=\overline{\mathrm{0,}N-1},$следует:

$\begin{array}{c}\xi S\left(x_{k+1}\right)=\frac{Y_{k+1}-Y_k}{h}-\frac{h_{}}{6}\left[m_k+2m_{k+1}\right],x\in \left[x_k,x_{k+1}\right],k=\overline{\mathrm{0,}N-1},\\ \xi S\left(0\right)=\frac{Y_1-Y_0}{h}+\frac{h_{}}{6}\left[2m_0+m_1\right].\end{array}$ 

Другой вычислительный алгоритм основывается на сведение модельной задачи (12) к системе интегральных уравнений. Решая дифференциальное уравнение (12), получили:

${\theta }_y\left(t,x\right)=e^{-t}{\theta }_{y_0}\left(x\right)+\underset{0}{\overset{t}{\int }}{e}^{-\left(t-\zeta \right)}\left[y\left(\zeta ,x\right)-f^{\left(y\right)}\left(\zeta ,x\right)\right]d\zeta ,$                                   (18)

где ${\theta }_{y_0}\left(x\right)\text{=}{y}_0\left(x\right)-{\xi }^2{y}_0\left(x\right)$.

Добавляя к выражению ${\theta }_y\left(t,x\right)\text{=}y\left(t,x\right)-{\xi }^{2}y\left(t,x\right)$ условие $y=0\text{ (}x=\mathrm{0,}x\text{=}\mathrm{1,}t>0),$ пришли к задаче, зависящей от параметра t:

$\left\{\begin{array}{l}{\xi }^{2}y\left(t,x\right)-y\left(t,x\right)=-{\theta }_y\left(t,x\right),\\ y\left(t\mathrm{,0}\right)=y\left(t\mathrm{,1}\right)=0.\end{array}\right.$                                                                 (19)

Решение, которой представлено в виде интегрального уравнения (20):

$y\left(t,x\right)=-\underset{0}{\overset{1}{\int }}G\left(x,s\right){\theta }_y\left(t,s\right)ds,$                                                                  (20)

$G\left(x,s\right)=\left\{\begin{array}{l}\frac{\mathrm{sh}s\mathrm{sh}\left(x-1\right)}{\mathrm{sh}1}\mathrm{,0}<s\le x,\\ \frac{\mathrm{sh}x\mathrm{sh}\left(s-1\right)}{\mathrm{sh}1},x<s\le x.\end{array}\right.$               

$\xi y\left(t,x\right)=-\underset{0}{\overset{x}{\int }}ch\left(x-s\right){\theta }_{\phi }\left(t,s\right)ds-\frac{chx}{sh1}\underset{0}{\overset{1}{\int }}ch\left(s-1\right){\theta }_{\phi }\left(t,s\right)ds$                         (21)

Проведя дискретизацию по переменной t в (18), (20), (21) и считая, что

${\left({\theta }_y\right)}_i^n\text{=}{\theta }_y\left(n\Delta ,x_i\right),y_i^n=y\left(n\Delta ,x_i\right),$ получили систему:

$\left\{\begin{array}{l}{y}_i^n=-\underset{0}{\overset{1}{\int }}G\left(x_i,s\right){\theta }_y\left(\left(n-1\right)\Delta ,s\right)ds,\\ \xi y\left(n\Delta ,x_i\right)=-\underset{0}{\overset{x_i}{\int }}ch\left(x_i-s\right){\theta }_{\phi }\left(\left(n-1\right)\Delta ,s\right)ds-\frac{ch{x}_i}{sh1}\underset{0}{\overset{1}{\int }}ch\left(s-1\right){\theta }_{\phi }\left(\left(n-1\right)\Delta ,s\right)ds,\\ {\left({\theta }_y\right)}_i^n=e^{-n\Delta }\left(\theta \left(\mathrm{0,}{x}_i\right)+\underset{0}{\overset{n\Delta }{\int }}{e}^{\zeta }\left[y\left(\zeta ,x_i\right)-f^{\left(y\right)}\left(\zeta ,x_i\right)\right]d\zeta \right),\end{array}\right.$       (22)

${\theta }_y\left(\mathrm{0,}{x}_{}\right)\text{=}{y}_0\left(x_i\right)-{\xi }^2{y}_0\left(x_i\right), n=\mathrm{1,2,....}$,

Расчет интегралов в (22) осуществлялся с помощью метода трапеции [9] с переменным шагом $h_i=x_{i+1}-x_i,i=\overline{\mathrm{0,}N-1}$ на сетке по x.

Переходим к рассмотрению последнего вычислительного алгоритма, который связан с использованием схемы предиктор−корректор. Для этого, в уравнении (12) и интегральных соотношениях провели дискретизацию по переменной t. В результате уравнение (12) записалось так:

$\frac{{\theta }_y^{*}-{\theta }_y^n}{\Delta /2}\text{+}{\theta }_y^{*}\text{=}{\left({\tilde{f}}^{\left(y\right)}\right)}^n,\frac{{\theta }_y^{n+1}-{\theta }_y^n}{\Delta }\text{+}{\theta }_y^{n+1}\text{=}{\left({\tilde{f}}^{\left(y\right)}\right)}^{*}.$

Из полученных уравнений выражаем ${\theta }_y^{*}, {\theta }_y^{n+1}$ :

${\theta }_y^{*}=\frac{{\theta }_y^n+\frac{\Delta }{2}{\left({\tilde{f}}^{\left(y\right)}\right)}^n}{1+\frac{\Delta }{2}},{\theta }_y^{n+1}=\frac{\Delta {\left({\tilde{f}}^{\left(y\right)}\right)}^{*}+{\theta }_y^n}{1+\Delta }.$

Используя интегральные уравнения (20) − (21), удалось перейти к следующей системе последовательных выражений:

$\left\{\begin{array}{l}{\left({\tilde{f}}^{\left(y\right)}\right)}_i^n=y_i^n-{\tilde{f}}^{\left(y\right)}\left(n\Delta ,x_i\right),{\left({\theta }_y^{*}\right)}_i=\frac{{\left({\theta }_y^{}\right)}_i^n+\frac{\Delta }{2}{\left({\tilde{f}}^{\left(y\right)}\right)}_i^n}{1+\frac{\Delta }{2}},\\ y_i^{*}=-\underset{0}{\overset{1}{\int }}G\left(x_i,s\right){\theta }_y^{*}\left(s\right)ds,\\ {\left({\tilde{f}}^{\left(y\right)}\right)}_i^{*}=y_i^{*}-{\tilde{f}}^{\left(y\right)}\left(n\Delta ,x_i\right),{\left({\theta }_y^{}\right)}_i^{n+1}=\frac{\Delta {\left({\tilde{f}}^{\left(y\right)}\right)}_i^{*}+{\left({\theta }_y^{}\right)}_i^n}{1+\Delta },\\ y_i^{n+1}=-\underset{0}{\overset{1}{\int }}G\left(x_i,s\right){\theta }_y^{}\left(\left(n+1\right)\Delta ,s\right)ds,\\ \xi y\left(\left(n+1\right)\Delta ,x_i\right)=-\underset{0}{\overset{x_i}{\int }}ch\left(x_i-s\right){\theta }_{\phi }\left(\left(n+1\right)\Delta ,s\right)ds-\\ -\frac{ch{x}_i}{sh1}\underset{0}{\overset{1}{\int }}ch\left(s-1\right){\theta }_{\phi }\left(\left(n+1\right)\Delta ,s\right)ds.\end{array}\right.$                  (23)

Расчет интегралов в (23) осуществлялся с помощью метода трапеции или метода Симпсона [9].

Таким образом, чтобы найти решение на следующем временном слое, необходимо задать начальные данные для E, R, φ, и последовательно осуществить действия:

1. Рассчитать правую часть ${\mathcal{F}}^{\left(\phi \right)}(\chi \left(n\Delta ,x\right),\rho )$ для первого уравнения в (9), где значение $\chi$ берется с n-го (предыдущего) слоя. Найти решение для электрического потенциала ${\phi }^{n+1}=\phi \left(\left(n+1\right)\Delta ,x\right)$ и вычислить производную $\xi {\phi }^{n+1}$.
2. Рассчитать правую часть для второго уравнения в (9) ${\mathcal{F}}^{\left(\sigma \right)}\left(\xi \sigma \left(n\Delta ,x\right),\xi \chi \left(n\Delta ,x\right),\xi \phi \left(\left(n+1\right)\Delta ,x\right)+B,{\sigma }^n\right)$, где σ, χ берутся с n-го слоя, а φ – с (n +1)-го слоя. Найти решение ${\sigma }^{n+1}=\sigma \left(\left(n+1\right)\Delta ,x\right)$ и вычислить производную $\xi {\sigma }^{n+1}$.
3. Рассчитать правую часть для третьего уравнения в (9) ${\mathcal{F}}^{\left(\chi \right)}\left(\xi \sigma \left(\left(n+1\right)\Delta ,x\right),\xi \chi \left(n\Delta ,x\right),\xi \phi \left(\left(n+1\right)\Delta ,x\right)+B,{\sigma }^{n+1},{\chi }^n,\rho \right),$ где значение χ берется с n-го слоя, а φ, σ – с (n +1)-го слоя. Найти решение ${\chi }^{n+1}=\chi \left(\left(n+1\right)\Delta ,x\right)$ и вычислить производную $\xi {\chi }^{n+1}$.

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

Для реализации и тестирования описанных вычислительных схем, разработана система «Численный анализ задачи о баллистическом диоде». Меню системы представлено на рисунке 2.

Рисунок 2 – Меню системы

В окне «Численное решение задачи о баллистическом диоде» (рисунок 3) необходимо задать физические и численные параметры, утвердить начальные данные или считать данные из файла, указать количество итераций, выбрать алгоритм решения задачи. Только после этого можно будет получить решение задачи и сохранить в файл. Более того, в системе возможно следить за поведением графика зависимости энергии от времени.

Расчеты показали, что разработанные вычислительные алгоритмы имеют разные скорости сходимости. Алгоритм, использующий кубические сплайны, имеет более высокую скорость, чем алгоритм, использующий схему предиктор-корректор.

Хорошие результаты при σ = 0,004 можно получить с помощью системы интегральных уравнений (рисунок 4). При небольших σ лучше использовать другие алгоритмы.

Проверку правильности работы разработанных алгоритмов осуществили с помощью метода ортогональной прогонки. Данный метод позволяет проводить расчеты с заданной точностью. На рисунке 5 представлены результаты метода ортогональной прогонки и разработанных алгоритмов. Из рисунка видно, что результаты практически совпадают, т.е. графики, построенные для плотности, потенциала и энергии электронов разными алгоритмами совпали. Есть небольшие расхождения только для скорости электронов.

Рисунок 3 – Окно «Численное решение задачи о баллистическом диоде»

Рисунок 4 – V = 1.5 Volt, σ = 0.004, L = 3 × 10^-7 m, N =100, Δ₁ = ⅟₁₂ , ε₁ = 10^-6.

Рисунок 5 – Расчеты четырех методов:

V = 1 Volt, σ = 0.2, L = 6 × 10^-7 m, N =200, Δ₁ = ⅟₁₂ , ε₁ = 10^-6.

Заключение и выводы

С помощью разработанной системы осуществлено огромное число расчетов при разных параметрах. Расчеты показали, что разработанные вычислительные алгоритмы дают физически правдоподобные результаты, которые совпадают с результатами других авторов. Материалы данной работы могут быть интересны специалистам, занимающихся моделированием различных физических процессов в полупроводниковых приборах.

About the authors

Alesya S. Shevchenko

Author for correspondence.
Email: ibragimova.a.s@mail.ru

Candidate of Physical and Mathematical Sciences, Associate Professor, Rubtsovsk Industrial Institute (branch) 

References

Supplementary files