半导体器件方程混合有限元三步两层网格法

Three-Step Two-Grid Algorithm Based on Mixed Finite Element Method for Semiconductor Device Equations

摘要: 半导体器件数值模拟一直是科学家们关注的重要领域。本文研究半导体器件的非线性方程组，对电子位势方程和两个浓度方程均采用混合有限元法进行离散，并构造了一种有效求解的三步两层网格算法。混合有限元法离散方程可以同时给出未知函数和未知函数通量同等阶数的误差逼近，也具有局部守恒性。三步两层网格算法不仅保持了数值解的可靠性和收敛阶，还大大缩短了计算时间。所以本研究具有重要的理论和实际意义，可为半导体器件的设计优化提供有力支持。

Abstract: Numerical simulation of semiconductor devices has always been an important area of interest for scientists. In this paper, the nonlinear equations of semiconductor device are studied. We use a mixed finite element method to discretize the existing electron potential equation and two concen-tration equations, and construct an effective three-step two-grid algorithm. The mixed finite ele-ment method for discretizing equations can simultaneously provide error approximations of the unknown function itself and its flux of the same order, and also has local conservation. The three-step two-grid algorithm not only maintains the reliability and convergence order of the nu-merical solution, but also greatly reduces computational time. The research has great theoretical and practical significance, providing strong support for the design and optimization of semiconduc-tor device.

1. 引言

用数值模拟技术研究半导体器件一直是科学家们关注的重要领域，因为半导体器件是各种集成电路和电子产品、光电探测器和发光二极管等光电器件、太阳能电池等绿色能源元件必不可少的重要组成部分，需求量巨大。然而，目前研究者遇到的难题是，不管用标准有限元法还是混合有限元法去离散半导体器件方程，均会得到超级大的代数方程组，直接影响求解时间，耗时漫长。

本文基于混合有限元法的离散，研究如何构造三步两层网格算法。混合有限元法的使用不仅可以同时给出方程组未知函数本身及其通量同等阶数的误差逼近，还具有局部守恒性。而两层网格算法是目前求解非线性方程组高效的算法之一，使用此算法既可以保持数值解的可靠性和收敛阶，又能降低计算时间。这对于推进数值模拟理论分析的发展，具有重大的理论和实际意义。

2. 半导体器件问题的数学模型

本文研究的问题是二维半导体器件中的电子位势和电子、空穴浓度的变化规律。该问题可以用一个包含三个偏微分方程的数学模型来描述，其中一个方程是椭圆型的电子位势方程，另外两个方程是抛物型的电子、空穴浓度方程。具体的方程形式如下 [1] ：

$\{\begin{array}{c}-\Delta \psi =\nabla \cdot u=\alpha \left(p-e+N\left(x\right)\right),\\ \frac{\partial e}{\partial t}=\nabla \cdot \left[D_e\left(x\right)\nabla e-{\mu }_e\left(x\right)e\nabla \psi \right]-R_{\text{1}}\left(e,p\right),\\ \frac{\partial p}{\partial t}=\nabla \cdot \left[D_p\left(x\right)\nabla p+{\mu }_p\left(x\right)p\nabla \psi \right]-R_{\text{2}}\left(e,p\right).\end{array}$ (1)

未知函数是电子位势 $\psi$ 以及电子浓度和空穴浓度e和p，方程的解具有一定的正则性。 $u=-\nabla \psi$ 是电场强度， $\alpha$ 是电子电荷与介电常数的比值，N(x)是净掺杂浓度， $D_e\left(x\right)$ 、 $D_p\left(x\right)$ 分别是电子扩散系数和空穴扩散系数， ${\mu }_e\left(x\right)$ 、 ${\mu }_p\left(x\right)$ 分别是电子移动率和空穴移动率， $R_i\left(e,p\right)$ (i =1,2)是关于e、p的产生复合率，这里的下标表示两个浓度方程的产生复合率可以不一样，比文献 [2] 中的研究更适用实际问题。其中 $\alpha$ 、D_e(x)、D_p(x)、𝜇_e(𝑥)、𝜇_p(𝑥)均有大于零的上、下界，此条件有利于估计数值解的误差范围，且 $x\in \Omega \subset R^{\text{2}}$ ， $t\in J=\left(0,T\right]$ 。我们考虑第一类边界条件： $\psi \left(x,t\right)=e\left(x,t\right)=p\left(x,t\right)=0$ ， $\left(x,t\right)\in \partial \Omega \times J$ ，第二、三类边界条件有待后期的研究。相应的初始条件为： $e\left(x,\text{0}\right)=e_{\text{0}}\left(x\right)$ ， $p\left(x,\text{0}\right)=p_{\text{0}}\left(x\right)$ 。

二维半导体器件模型描述了半导体器件中电子位势和载流子浓度随时间和空间的变化规律。为了求解上述方程组，本文采用混合有限元法离散电子位势方程和电子和空穴的两个浓度方程。这是一种常用的数值求解方法，其他常用的方法还包括差分法、有限体积法等。本文将在下一节给出混合有限元方法在求解半导体问题中的应用，并讨论该方法的优点。

3. 混合有限元法全离散格式

首先在 $L^2\left(\Omega \right)$ 空间的基础上定义空间 $W=\left\{w\in L^2\left(\Omega \right),\left(w,1\right)=0\right\}$ 和 $V=\left\{v\in {\left(L^2\left(\Omega \right)\right)}^2,\nabla \cdot v\in L^2\left(\Omega \right)\right\}$ 。令 $z_e=D_e\nabla e,z_p=D_p\nabla p$ ，则两者均是与未知量相关的变量。方程组(1)的弱形式为在空间 $W\times V\times L^2\left(\Omega \right)\times V\times L^2\left(\Omega \right)\times V$ 上找满足方程组

$\{\begin{array}{c}\left(u,v\right)-\left(\nabla \cdot v,\psi \right)=0,\forall v\in V,\\ \left(\nabla \cdot u,w\right)=\alpha \left(p-e+N,w\right),\forall w\in W,\\ \begin{array}{l}\left(\frac{\partial e}{\partial t},\phi \right)-\left(\nabla \cdot z_e,\phi \right)-\left(eu\cdot \nabla {\mu }_e,\phi \right)-\left({\mu }_e\nabla e\cdot u,\phi \right)\\ -\alpha \left({\mu }_{e}e\left(p-e+N\right),\phi \right)=-\left(R_{\text{1}}\left(e,p\right),\phi \right),\forall \phi \in L^2\left(\Omega \right),\end{array}\\ \left(D_e^{-1}{z}_e,\chi \right)+\left(e,\nabla \cdot \chi \right)=0,\forall \chi \in V,\\ \begin{array}{l}\left(\frac{\partial p}{\partial t},\phi \right)-\left(\nabla \cdot z_p,\phi \right)+\left(pu\cdot \nabla {\mu }_p,\phi \right)+\left({\mu }_p\nabla p\cdot u,\phi \right)\\ +\alpha \left({\mu }_{p}p\left(p-e+N\right),\phi \right)=-\left(R_{\text{2}}\left(e,p\right),\phi \right),\forall \phi \in L^2\left(\Omega \right),\end{array}\\ \left(D_p^{-1}{z}_p,\chi \right)+\left(p,\nabla \cdot \chi \right)=0,\forall \chi \in V.\end{array}$ (2)

的弱解 $\left(\psi ,u,e,z_e,p,z_p\right)$ 。

对电子位势方程采用混合有限元方法，选择网格尺寸为 $h_{\psi }$ 的有限元剖分 ${\Gamma }_{h\psi }$ ，取包含于空间 $W\times V$ 的 $W_h\times V_{h\psi }$ 作为混合有限元子空间，其阶为k。对两个浓度方程同样采用混合有限元方法，选择网格尺寸为 $h_c$ 的有限元剖分 ${\Gamma }_{hc}$ ，取包含于空间 $L^2\left(\Omega \right)\times V$ 的 $M_h\times V_{hc}$ 作为混合有限元子空间，其阶为l。对时间进行正常的剖分，则方程组(2)可转化为在空间 $W_h\times V_{h\phi }\times M_h\times V_{hc}\times M_h\times V_{hc}$ 上找满足方程组

$\{\begin{array}{c}\left(u_h^n,v\right)-\left(\nabla \cdot v,{\psi }_h^n\right)=0,\forall v\in V_{h\psi },\\ \left(\nabla \cdot u_h^n,w\right)=\alpha \left(p_h^n-e_h^n+N,w\right),\forall w\in W_h,\\ \begin{array}{l}\left(\frac{e_h^n-e_h^{n-1}}{\Delta t},\phi \right)-\left(\nabla \cdot z_{eh}^n,\phi \right)-\left(e_h^n{u}_h^n\cdot \nabla {\mu }_e,\phi \right)-\left({\mu }_e\nabla e_h^n\cdot u_h^n,\phi \right)\\ -\alpha \left({\mu }_e{e}_h^n\left(p_h^n-e_h^n+N\right),\phi \right)=-\left(R_1\left(e_h^n,p_h^n\right),\phi \right),\forall \phi \in M_h,\end{array}\\ \left(D_e^{-1}{z}_{eh}^n,\chi \right)+\left(e_h^n,\nabla \cdot \chi \right)=0,\forall \chi \in V_{hc}\\ \begin{array}{l}\left(\frac{p_h^n-p_h^{n-1}}{\Delta t},\phi \right)-\left(\nabla \cdot z_{ph}^n,\phi \right)+\left(p_h^n{u}_h^n\cdot \nabla {\mu }_p,\phi \right)+\left({\mu }_p\nabla p_h^n\cdot u_h^n,\phi \right)\\ +\alpha \left({\mu }_p{p}_h^n\left(p_h^n-e_h^n+N\right),\phi \right)=-\left(R_2\left(e_h^n,p_h^n\right),\phi \right),\forall \phi \in M_h.\end{array}\\ \left(D_p^{-1}{z}_{ph}^n,\chi \right)+\left(p_h^n,\nabla \cdot \chi \right)=0,\forall \chi \in V_{hc}.\end{array}$ (3)

的解 $\left({\psi }_h^n,u_h^n,e_h^n,z_{eh}^n,p_h^n,z_{ph}^n\right)$ 。其中浓度方程的初始值 $e_h^0,p_h^0$ 由它们的投影决定。

定理1 已知 $\left({\psi }^n,u^n,e^n,z_e^n,p^n,z_p^n\right)$ 是方程组(2)的弱解， $\left({\psi }_h^n,u_h^n,e_h^n,z_{eh}^n,p_h^n,z_{ph}^n\right)$ 是方程组(3)的有限元数值解，则对 $1\le n\le N$ 和 $l,k\ge 1$ ，有收敛结果

${\Vert {\psi }^n-{\psi }_h^n\Vert }_W+{\Vert u^n-u_h^n\Vert }_V\le C\left(h_c^{l+1}+h_{\psi }^{k+1}+\Delta t\right),$

${\Vert s^n-s_h^n\Vert }_{0,\text{4}}+{\Vert z_s^n-z_{sh}^n\Vert }_{0,\text{4}}\le C\left(h_c^{l+1}+h_{\psi }^{k+1}+\Delta t\right),s=e,p.$

定理1 的证明可参考文献 [2] 中5.2混合有限元—特征混合有限元法数值解的误差估计。

根据方程组(3)和定理1可知，采用混合有限元法离散方程，最后可以同时逼近两个不同的未知量，其中一个未知量与另一个未知变量的导数有物理上的关联。例如，电场强度 $u=-\nabla \psi$ 在实际应用中是一种非常有用的物理量。如果使用标准有限元法，需要先求解电子位势 $\psi$ ，然后进行微分运算才能得到电场强度的数值解。但是，采用混合有限元法可以同时解出电子位势 $\psi$ 和电场强度u的数值解，从而保证在相同的计算工作量下获得更高的精度。而且电子位势 $\psi$ 和电场强度u的误差收敛阶是一致的。此外，混合有限元法还具有局部守恒性。因此，采用混合有限元法来求解二维半导体器件模型方程组是一种非常可行的方法 [2] 。

4. 三步两层网格算法的分析与构造

由混合有限元法离散半导体器件方程组，在计算中会生成一个非常庞大的系数矩阵，导致半导体器件数值模拟计算是超大规模的问题，求解所需时间漫长，所以急需构造一个既能保证数值解的可靠性和收敛阶，又能降低计算时间的有效求解算法。我们选用的是两层网格算法。两层网格法最初是由许进超教授提出用于解决非线性椭圆型问题 [3] ，现在已成为高效求解各类复杂问题，如非线性方程组问题的重要方法之一。该方法的基本思想是设置两种不同尺寸的网格步长，例如较大较粗的网格步长H和较小较细的网格步长h，先后在两个不同网格步长的子空间上求解方程。即先在粗网格上求解原有给定的方程(一般都是较复杂的，求解特别耗时的方程组)，得到一组数值解，因为此时的网格步长较大较粗，所以求解过程所耗时间不多；然后将原有方程线性化，以粗网格上得到的数值解为初始解，在细网格上求解已线性化的方程(相比非线性方程组，此时的较简单，求解过程所耗时间较少)，以此达到降低计算时间的目的 [4] [5] 。

我们针对半导体器件方程组，构造三步两层网格算法。设 ${\Gamma }_H$ 是网格步长为H的粗网格， ${\Gamma }_h$ 是网格步长为h的细网格，且 $h\ll H<1$ 。定义空间 $W_H\times V_{H\phi }\times M_H\times V_{Hc}\times M_H\times V_{Hc}$ 是包含于 $W_h\times V_{h\phi }\times M_h\times V_{hc}\times M_h\times V_{hc}$ 上的粗空间。首先在粗网格上求解原有的非线性方程组，得到粗空间上初始逼近解 $\left({\psi }_H^n,u_H^n,e_H^n,z_{eH}^n,p_H^n,z_{pH}^n\right)$ 。然后在细网格上求解基于牛顿迭代线性化的方程组，得到两步两层网格解 $\left({\Psi }_h^n,U_h^n,E_h^n,Z_{eh}^n,P_h^n,Z_{ph}^n\right)$ 。这里是通过对乘积未知项的分解，获得线性化的结果，例如 $E_h^n{U}_h^n\approx e_H^n{U}_h^n+\left(E_h^n-e_H^n\right)u_H^n$ 。最后我们在细网格上再次求解线性化的方程组，得到校正的三步两层网格解 $\left({\tilde{\psi }}_h^n,{\tilde{u}}_h^n,{\tilde{e}}_h^n,{\tilde{z}}_{eh}^n,{\tilde{p}}_h^n,{\tilde{z}}_{ph}^n\right)$ 。具体步骤如下：

步骤一：在 ${\Gamma }_H$ 上求解下述原始的非线性方程组的解 $\left({\psi }_H^n,u_H^n,e_H^n,z_{eH}^n,p_H^n,z_{pH}^n\right)$

$\{\begin{array}{c}\left(u_H^n,v\right)-\left(\nabla \cdot v,{\psi }_H^n\right)=0,\forall v\in V_{H\psi },\\ \left(\nabla \cdot u_H^n,w\right)=\alpha \left(p_H^n-e_H^n+N,w\right),\forall w\in W_H,\\ \begin{array}{l}\left(\frac{e_H^n-e_H^{n-1}}{\Delta t},\phi \right)-\left(\nabla \cdot z_{eH}^n,\phi \right)-\left(e_H^n{u}_H^n\cdot \nabla {\mu }_e,\phi \right)-\left({\mu }_e\nabla e_H^n\cdot u_H^n,\phi \right)\\ -\alpha \left({\mu }_e{e}_H^n\left(p_H^n-e_H^n+N\right),\phi \right)=-\left(R_{\text{1}}\left(e_H^n,p_H^n\right),\phi \right),\forall \phi \in M_H,\end{array}\\ \left(D_e^{-1}{z}_{eH}^n,\chi \right)+\left(e_H^n,\nabla \cdot \chi \right)=0,\forall \chi \in V_{Hc}\\ \begin{array}{l}\left(\frac{p_H^n-p_H^{n-1}}{\Delta t},\phi \right)-\left(\nabla \cdot z_{pH}^n,\phi \right)+\left(p_H^n{u}_H^n\cdot \nabla {\mu }_p,\phi \right)+\left({\mu }_p\nabla p_H^n\cdot u_H^n,\phi \right)\\ +\alpha \left({\mu }_p{p}_H^n\left(p_H^n-e_H^n+N\right),\phi \right)=-\left(R_{\text{2}}\left(e_H^n,p_H^n\right),\phi \right),\forall \phi \in M_H.\end{array}\\ \left(D_p^{-1}{z}_{pH}^n,\chi \right)+\left(p_H^n,\nabla \cdot \chi \right)=0,\forall \chi \in V_{Hc}.\end{array}$ (4)

步骤二：在 ${\Gamma }_h$ 上求解下述已经线性化的方程组的解 $\left({\Psi }_h^n,U_h^n,E_h^n,Z_{eh}^n,P_h^n,Z_{ph}^n\right)$

$\{\begin{array}{c}\left(U_h^n,v\right)-\left(\nabla \cdot v,{\Psi }_h^n\right)=0,\forall v\in V_{h\phi },\\ \left(\nabla \cdot U_h^n,w\right)=\alpha \left(P_h^n-E_h^n+N,w\right),\forall w\in W_h,\\ \begin{array}{l}\left(\frac{E_h^n-E_h^{n-1}}{\Delta t},\phi \right)-\left(\nabla \cdot Z_{eh}^n,\phi \right)-\left(\left[e_H^n{U}_h^n+\left(E_h^n-e_H^n\right)u_H^n\right]\cdot \nabla {\mu }_e,\phi \right)\\ -\left({\mu }_e\left[\nabla e_H^n\cdot U_h^n+\nabla \left(E_h^n-e_H^n\right)\cdot u_H^n\right],\phi \right)-\alpha ({\mu }_e[e_H^n\left(P_h^n-E_h^n+N\right)\\ +\left(E_h^n-e_H^n\right)\left(p_H^n-e_H^n+N\right)],\phi )=-\left({\text{R}}_{\text{1}}\left(E_h^n,P_h^n\right),\phi \right),\forall \phi \in M_h,\end{array}\\ \left(D_e^{-1}{Z}_{eh}^n,\chi \right)+\left(E_h^n,\nabla \cdot \chi \right)=0,\forall \chi \in V_{hc}\\ \begin{array}{l}\left(\frac{P_h^n-P_h^{n-1}}{\Delta t},\phi \right)-\left(\nabla \cdot Z_{ph}^n,\phi \right)+\left(\left[p_H^n{U}_h^n+\left(P_h^n-p_H^n\right)u_H^n\right]\cdot \nabla {\mu }_p,\phi \right)\\ +\left({\mu }_p\left[\nabla p_H^n\cdot U_h^n+\nabla \left(P_h^n-p_H^n\right)\cdot u_H^n\right],\phi \right)+\alpha ({\mu }_p[p_H^n\left(P_h^n-E_h^n+N\right)\\ +\left(P_h^n-p_H^n\right)\left(p_H^n-e_H^n+N\right)],\phi )=-\left({\text{R}}_{\text{2}}\left(E_h^n,P_h^n\right),\phi \right),\forall \phi \in M_h\end{array}\\ \left(D_p^{-1}{Z}_{ph}^n,\chi \right)+\left(P_h^n,\nabla \cdot \chi \right)=0,\forall \chi \in V_{hc}\end{array}$ (5)

其中对于产生复合率的处理如下：

${\text{R}}_{\text{1}}\left(E_h^n,P_h^n\right)=R_{\text{1}}\left(e_H^n,p_H^n\right)+R_{\text{1}e}\left(e_H^n,p_H^n\right)\left(E_h^n-e_H^n\right)+R_{\text{1}p}\left(e_H^n,p_H^n\right)\left(P_h^n-p_H^n\right)$ ，

上式中的下标e或p表示函数 $R_{\text{1}}\left(e^n,p^n\right)$ 对变量的偏导，类似可定义 ${\text{R}}_{\text{2}}\left(E_h^n,P_h^n\right)$ 。

步骤三：在 ${\Gamma }_h$ 上求解下述线性方程组的解 $\left({\tilde{\psi }}_h^n,{\tilde{u}}_h^n,{\tilde{e}}_h^n,{\tilde{z}}_{eh}^n,{\tilde{p}}_h^n,{\tilde{z}}_{ph}^n\right)$

$\{\begin{array}{c}\left({\tilde{u}}_h^n,v\right)-\left(\nabla \cdot v,{\tilde{\psi }}_h^n\right)=0,\forall v\in V_{h\phi },\\ \left(\nabla \cdot {\tilde{u}}_h^n,w\right)=\alpha \left({\tilde{p}}_h^n-{\tilde{e}}_h^n+N,w\right),\forall w\in W_h,\\ \begin{array}{l}\left(\frac{{\tilde{e}}_h^n-{\tilde{e}}_h^{n-1}}{\Delta t},\phi \right)-\left(\nabla \cdot {\tilde{z}}_{eh}^n,\phi \right)-\left(\left[E_h^n{\tilde{u}}_h^n+\left({\tilde{e}}_h^n-E_h^n\right)U_h^n\right]\cdot \nabla {\mu }_e,\phi \right)\\ -\left({\mu }_e\left[\nabla E_h^n\cdot {\tilde{u}}_h^n+\nabla \left({\tilde{e}}_h^n-E_h^n\right)\cdot U_h^n\right],\phi \right)-\alpha ({\mu }_e[E_h^n\left({\tilde{p}}_h^n-{\tilde{e}}_h^n+N\right)\\ +\left({\tilde{e}}_h^n-E_h^n\right)\left(P_h^n-E_h^n+N\right)],\phi )=-\left({\Re }_1\left({\tilde{e}}_h^n,{\tilde{p}}_h^n\right),\phi \right),\forall \phi \in M_h,\end{array}\\ \left(D_e^{-1}{\tilde{z}}_{eh}^n,\chi \right)+\left({\tilde{e}}_h^n,\nabla \cdot \chi \right)=0,\forall \chi \in V_{hc}\\ \begin{array}{l}\left(\frac{{\tilde{p}}_h^n-{\tilde{p}}_h^{n-1}}{\Delta t},\phi \right)-\left(\nabla \cdot {\tilde{z}}_{ph}^n,\phi \right)+\left(\left[P_h^n{\tilde{u}}_h^n+\left({\tilde{p}}_h^n-P_h^n\right)U_h^n\right]\cdot \nabla {\mu }_p,\phi \right)\\ +\left({\mu }_p\left[\nabla P_h^n\cdot {\tilde{u}}_h^n+\nabla \left({\tilde{p}}_h^n-P_h^n\right)\cdot U_h^n\right],\phi \right)+\alpha ({\mu }_p[P_h^n\left({\tilde{p}}_h^n-{\tilde{e}}_h^n+N\right)\\ +\left({\tilde{p}}_h^n-P_h^n\right)\left(P_h^n-E_h^n+N\right)],\phi )=-\left({\Re }_2\left({\tilde{e}}_h^n,{\tilde{p}}_h^n\right),\phi \right),\forall \phi \in M_h\end{array}\\ \left(D_p^{-1}{\tilde{z}}_{ph}^n,\chi \right)+\left({\tilde{p}}_h^n,\nabla \cdot \chi \right)=0,\forall \chi \in V_{hc}\end{array}$ (6)

其中 ${\Re }_{\text{1}}\left({\tilde{e}}_h^n,{\tilde{p}}_h^n\right)=R\left(E_h^n,P_h^n\right)+R_{\text{1}e}\left(E_h^n,P_h^n\right)\left({\tilde{e}}_h^n-E_h^n\right)+R_{\text{1}p}\left(E_h^n,P_h^n\right)\left({\tilde{p}}_h^n-P_h^n\right)$ 。同样，这里的下标e或p表示函数 $R_{\text{1}}\left(e^n,p^n\right)$ 对变量的偏导，类似可定义 ${\Re }_{\text{2}}\left({\tilde{e}}_h^n,{\tilde{p}}_h^n\right)$ 。

5. 结束语

本文研究的半导体器件数值模拟问题，在分析半导体物理性能、设计和优化超大规模集成电路、研究新型结构器件等多个方面均有参与，也是促进发展我国半导体芯片行业的一个重要环节。在本文中，对电子位势方程和电子及空穴的两个浓度方程均采用混合有限元法进行离散，可以同时给出未知函数本身及其通量同等阶数的误差逼近，而混合有限元也具有局部质量守恒这个非常好的属性，又没有特征有限元复杂的分析和计算过程。为了减少计算时间，本文提出了一种三步两层网格算法，先在粗网格上求解非线性方程组，然后在细网格上求解线性化的方程组，而不是直接在细网格上求解大规模的原始方程。因为粗网格空间维数比细网格空间的小很多，所以可大大减少计算时间。并且增加的第三步可以进行数值解的校正，提高收敛阶，参见文献 [2] 。不足之处在于未能进行数值实验以验证理论分析结果，这将是我们后期的研究工作。

基金项目

2021年度国家自然科学基金青年科学基金项目“半导体器件问题的几类高效有限元两层网格法”(12101224)，2021年湖南农业大学大学生创新训练项目(XCX2021002)。

参考文献