用于半导体连续性方程的流线迎风有限元方法及系统

摘要

本发明公开了一种用于半导体连续性方程的混合流线迎风有限元方法及系统，其方法包括如下步骤：S1、几何模型的空间离散化；S2、网格单元的构造及应用；S3、利用Nedelec边缘基空间将电流密度模型在网格单元相邻段的中心进行插值得到电流密度；S4、单元矩阵方程的构造；S5、通过对求解域中的所有网格单元进行遍历，得到系统矩阵方程。相较于FBSG方法，本发明对空间网格质量要求更低，从而使得收敛性更好；相较于SUPG方法，本发明在迎风函数的多维应用方面进行了优化；相较于FVFEM‑SG方法，本发明在定义人工扩散系数方面更具灵活性。

权利要求

1.一种用于半导体连续性方程的混合流线迎风有限元方法，其特征在于，包括如下步骤：

S1、几何模型的空间离散化：将物理模型简化为每个网格上的时域差分方程，而时域进行求解取决于器件的运行方式，在静态直流条件下，时间差设置为零；

S2、网格单元的构造及应用：所述网格单元为四面体形，由网格顶点、网格边缘以及网格面组成；网格单元的应用是将电子连续性方程在网格单元上进行积分；步骤S2中，网格单元的应用包括：

S2.1.将电子连续性方程在网格单元上进行加权积分：其中q为单位电荷，n为电子浓度，Jn为电子电流密度，Ωe为网格单元，Rn为净复合率，Nk为加权函数，dΩ为体积微元； 为电子浓度随时间的变化率；

S2.2.通过标准拉格朗日插值函数得出每个网格内的电子密度：m∈Ωe∪ΓN是指顶点vm位于计算域内或者纽曼边界上，m∈ΓD是指顶点vm属于狄利克雷边界，nd,m(t)表示相应的与时间相关的狄利克雷边界值，Nm为插值展开函数，nm为节点上的电子浓度；

S2.3.沿网格边缘eij的SU电流密度模型表示为：其中 是沿边缘eij的电场，μn,ij是边缘电子迁移率，Dn,ij是边缘扩散系数， 是辅助人工扩散系数，其在算法性能方面起着关键作用；另外辅助人工扩散系数 满足关系式其中lij为边缘eij的长度，α为迎风系数， 为佩克莱特数，其中 为电势，μn以及Dn分别为载流子迁移率、载流子扩散率；

S3、利用Nedelec边缘基空间将式(3)中的SU电流密度模型在 相邻段的中心进行插值，得到

其中 为电流密度， 为矢量基函数S4、单元矩阵方程的构造，具体包括：S4.1.通过将公式(2)、(5)代入公式(1)并在整个Ωe区域进行积分，得到：其中h为纽曼边界条件；

S4.2.尽管公式(6)左边部分的第二项是在网格边缘所进行的操作，但是节点的未知数仍被计算，通过将边缘电流密度函数公式(3)代入公式(6)，因此公式(6)又写为：其中当定点为i和j时，σr分别为‑1和1；

S4.3.通过遍历网格Ωe中所有顶点，得到公式(7)的矩阵形式：[Ke]是单位“迎风”刚度矩阵，其单元由公式(9)计算出[Me]是阻尼矩阵，其元素表示为：fe是列向量，其元素表示为：在应用时域后向差分之后，公式(8)离散为：([Me]‑[Ke]){n}t+Δt＝[Me]{n}t+{fe}   (12)其中Δt是演化时间的步长；

S5、通过对求解域中的所有网格单元进行遍历，得到系统矩阵方程。

2.根据权利要求1所述用于半导体连续性方程的混合流线迎风有限元方法，其特征在于，步骤S5中，系统矩阵方程的构造具体如下：通过对求解域中的所有网格单元进行遍历，得出系统矩阵方程：([M]‑[K]){n}t+Δt＝[M]{n}t+{f}    (13)其中[K]是系统“迎风”刚度矩阵，[M]是系统阻尼矩阵；在稳态情况下，与时间有关的项均应该设置为0，得到：

‑[K]{n}＝{f}   (14)公式(13)以及公式(14)两个矩阵方程分别用于求解瞬态情况和稳态情况。

3.一种用于半导体连续性方程的混合流线迎风有限元系统，其特征在于，包括如下模块：

几何模型的空间离散化模块：将物理模型简化为每个网格上的时域差分方程，而时域进行求解取决于器件的运行方式，在静态直流条件下，时间差设置为零；

网格单元的构造及应用模块：所述网格单元为四面体形，由网格顶点、网格边缘以及网格面组成；网格单元的应用是将电子连续性方程在网格单元上进行积分；网格单元的应用包括：

将电子连续性方程在网格单元上进行加权积分：其中q为单位电荷，n为电子浓度，Jn为电子电流密度，Ωe为网格单元，Rn为净复合率，Nk为加权函数，dΩ为体积微元； 为电子浓度随时间的变化率；

通过标准拉格朗日插值函数得出每个网格内的电子密度：m∈Ωe∪ΓN是指顶点vm位于计算域内或者纽曼边界上，m∈ΓD是指顶点vm属于狄利克雷边界，nd,m(t)表示相应的与时间相关的狄利克雷边界值，Nm为插值展开函数，nm为节点上的电子浓度；

格边缘eij的SU电流密度模型表示为：其中 是沿边缘eij的电场，μn,ij是边缘电子迁移率，Dn,ij是边缘扩散系数， 是辅助人工扩散系数，其在算法性能方面起着关键作用；另外辅助人工扩散系数 满足关系式其中lij为边缘eij的长度，α为迎风系数， 为佩克莱特数，其中 为电势，μn以及Dn分别为载流子迁移率、载流子扩散率；

电子电流密度计算模块：利用Nedelec边缘基空间将式(3)中的SU电流密度模型在相邻段的中心进行插值，得到

其中 为电流密度， 为矢量基函数；

单元矩阵方程构造模块：构造单元矩阵方程；单元矩阵方程的构造模块具体如下：通过将公式(2)、(5)代入公式(1)并在整个Ωe区域进行积分，得到：其中h为纽曼边界条件；

尽管公式(6)左边部分的第二项是在网格边缘所进行的操作，但是节点的未知数仍被计算，通过将边缘电流密度函数公式(3)代入公式(6)，因此公式(6)又写为：其中当定点为i和j时，σr分别为‑1和1；

通过遍历网格Ωe中所有顶点，得到公式(7)的矩阵形式：[Ke]是单位“迎风”刚度矩阵，其单元由公式(9)计算出[Me]是阻尼矩阵，其元素表示为：fe是列向量，其元素表示为：在应用时域后向差分之后，公式(8)离散为：([Me]‑[Ke]){n}t+Δt＝[Me]{n}t+{fe}  (12)其中Δt是演化时间的步长；

系统矩阵方程构造模块：通过对求解域中的所有网格单元进行遍历，得到系统矩阵方程。

4.根据权利要求3所述用于半导体连续性方程的混合流线迎风有限元系统，其特征在于，所述系统矩阵方程的构造模块具体如下：通过对求解域中的所有网格单元进行遍历，得出系统矩阵方程：([M]‑[K]){n}t+Δt＝[M]{n}t+{f}    (13)其中[K]是系统“迎风”刚度矩阵，[M]是系统阻尼矩阵；在稳态情况下，与时间有关的项均应该设置为0，得到：

‑[K]{n}＝{f}  (14)公式(13)以及公式(14)两个矩阵方程分别用于求解瞬态情况和稳态情况。

说明书

用于半导体连续性方程的流线迎风有限元方法及系统

技术领域

[0001] 本发明应用于半导体仿真技术领域，特别是半导体器件中以对流为主导的非线性载流子输运方程的数值离散化以及混合数值方法的构造和分析。

背景技术

[0002] 数值仿真技术作为一种实用的计算机辅助设计工具，其使得数字以及模拟应用半导体器件的设计以及优化更加便利。在固体电子器件的数学物理模型中，漂移扩散(Drift‑
Diffusion,DD)模型是最常用的。DD模型是基于电子和空穴电流连续性方程的漂移扩散图，
并与泊松方程耦合。然而，电子和空穴电流连续性方程的混合性质使得其具有很强的非线
性，这就使得其空间离散变得尤为重要。首先，研究人员试图使用基于传统有限差分方法
(Finite Differential Method,FDM)或有限元方法(Finite Element Method,FEM)的载流
子电流连续性方程。然而，在数值近似解中经常会发生非物理振荡，特别是在载流子密度梯
度很大的区域。经过数十年来广泛的研究，研究人员发现该稳定性可以通过对传统数值算
法增加额外的人工扩散系数来实现。然而，起辅助作用的人工扩散系数不仅要稳定数值算
法，而且不能降低数值求解精度。因此，找到一个最优人工扩散系数的有效结构是至关重要
的。目前，FB‑SG(Finite Box‑Scharfetter‑Gummel)算法在许多商业仿真软件中被广泛应
用，然而FB‑SG通常需要高网格质量的空间离散化。此外，SUPG(Streamline Upwind Petrov 
Galerkin)也是另一种有效方法，与SG模型相比，一维SU模型在人工扩散系数的定义方面更
具有灵活性。然而，其在二维以及三维应用中寻找“最佳”人工扩散系数通常会有很大的挑
战性。

发明内容

[0003] 针对上述问题，本发明提供了一种用于半导体连续性方程的混合流线迎风有限元方法(SU‑FEM)及系统，本发明方法相较于有限体积有限元SG法(FVFEM‑SG)，其在定义人工
扩散系数方面更具灵活性；相较于SUPG法，本发明方法在多维应用时可靠性更高。

[0004] 本发明采用如下技术方案：

[0005] 一种用于半导体连续性方程的混合流线迎风有限元方法，其步骤包括：

[0006] S1、几何模型的空间离散化；空间离散过程是将物理模型简化为每个网格上的时域差分方程，而时域求解取决于器件的运行方式，在静态直流条件下，可以直接将时差设置
为零；

[0007] S2、网格单元的构造及应用；所述网格单元为四面体形，由网格顶点、网格边缘以及网格面组成；网格单元的应用是将电子连续性方程在网格单元上进行积分；

[0008] S3、利用Nedelec边缘基空间将电流密度模型在网格单元相邻段的中心进行插值得到电流密度；

[0009] S4、单元矩阵方程的构造；

[0010] S5、通过对求解域中的所有网格单元进行遍历，得到系统矩阵方程。

[0011] 优选的，步骤2中，网格单元的具体应用如下：

[0012] S2.1.将电子连续性方程在网格单元上进行加权积分：

[0013]

[0014] 其中，q为单位电荷，n为电子浓度， 为电子电流密度，Ωe为网格单元，Rn为净复合率，Nk为加权函数，dΩ为体积微元；

[0015] S2.2.通过标准拉格朗日插值函数可以得出每个网格内的电子浓度：

[0016]

[0017] m∈Ωe∪ΓN是指顶点vm位于计算域内或者纽曼边界上，m∈ΓD是指顶点vm属于狄利克雷边界，nd,m(t)表示相应的与时间相关的狄利克雷边界值，Nm为插值展开函数，nm为节
点上的电子浓度；

[0018] S2.3.为了更清楚地描述混合数值算法的构造，下面将只阐述在网格Ωe上的积分。为处理连续性方程的非线性，引入了SU电流密度模型来处理电子和空穴电流连续性方
程。沿网格边缘eij的SU电流密度模型可以表示为：

[0019]

[0020] 其中 是沿边缘eij的电场，μn,ij是边缘电子迁移率，Dn,ij是边缘扩散系数，是辅助人工扩散系数，其在算法性能方面起着关键作用；另外辅助人工扩散系数 满足
关系式

[0021]

[0022] 其中lij为边缘eij的长度，α为迎风系数， 为佩克莱特数，其中 为电势，μn以及Dn分别为载流子迁移率、载流子扩散率。

[0023] 优选的，步骤3具体如下：

[0024] 混合流线迎风有限元方法通过利用Nedelec边缘基空间将式(3)中的SU电流密度模型在 相邻段的中心进行插值，因此可以得到

[0025]

[0026] 其中 为电流密度， 为矢量基函数。

[0027] 优选的，步骤4中，单元矩阵方程的构造过程具体如下：

[0028] S4.1.通过将公式(2)、(5)代入公式(1)并在整个Ωe区域进行积分，可以得到：

[0029]

[0030] 其中h为纽曼边界条件；

[0031] S4.2.尽管公式(6)的左边部分的第二项是在网格边缘所进行的操作，但是节点的未知数仍被计算，通过将边缘电流密度函数公式(3)代入公式(6)，因此公式(6)又可写为：

[0032]

[0033] 其中当定点为i和j时，σr分别为‑1和1；

[0034] S4.3.通过遍历网格Ωe中所有顶点，可以得到公式(7)的矩阵形式：

[0035]

[0036] [Ke]是单位“迎风”刚度矩阵，其单元由公式(9)计算出

[0037]

[0038] [Me]是阻尼矩阵，其元素可以表示为：

[0039]

[0040] fe是列向量，其元素可以表示为：

[0041]

[0042] 在应用时域后向差分之后，公式(8)可以离散为：

[0043] ([Me]‑[Ke]){n}t+Δt＝[Me]{n}t+{fe}  (12)

[0044] 其中Δt是演化时间的步长。

[0045] 优选的，步骤5中系统矩阵方程的形成具体如下：

[0046] S5.1.通过对求解域中的所有网格单元进行遍历，可以得出系统矩阵方程：

[0047] ([M]‑[K]){n}t+Δt＝[M]{n}t+{f}  (13)

[0048] 其中[K]是系统“迎风”刚度矩阵，[M]是系统阻尼矩阵。另外，在稳态情况下，与时间有关的项均应该设置为0，因此可以得到：

[0049] ‑[K]{n}＝{f}  (14)

[0050] 公式(13)以及公式(14)两个矩阵方程可以分别用于求解瞬态情况和稳态情况。

[0051] 本发明还公开了一种用于半导体连续性方程的混合流线迎风有限元系统，其包括如下模块：

[0052] 几何模型的空间离散化模块：将物理模型简化为每个网格上的时域差分方程，而时域进行求解取决于器件的运行方式，在静态直流条件下，时间差设置为零；

[0053] 网格单元的构造及应用模块：所述网格单元为四面体形，由网格顶点、网格边缘以及网格面组成；网格单元的应用是将电子连续性方程在网格单元上进行积分；

[0054] 电子电流密度计算模块：利用Nedelec边缘基空间将电流密度模型在网格单元相邻段的中心进行插值得到电流密度；

[0055] 单元矩阵方程构造模块：构造单元矩阵方程；

[0056] 系统矩阵方程构造模块：通过对求解域中的所有网格单元进行遍历，得到系统矩阵方程。

[0057] 优选的，所述网格单元的具体应用包括：

[0058] 将电子连续性方程在网格单元上进行加权积分：

[0059]

[0060] 其中q为单位电荷，n为电子浓度，Jn为电子电流密度，Ωe为网格单元，Rn为净复合率，Nk为加权函数，dΩ为体积微元； 为电子浓度随时间的变化率；

[0061] 通过标准拉格朗日插值函数得出每个网格内的电子密度：

[0062]

[0063] m∈Ωe∪ΓN是指顶点vm位于计算域内或者纽曼边界上，m∈ΓD是指顶点vm属于狄利克雷边界，nd,m(t)表示相应的与时间相关的狄利克雷边界值，Nm为插值展开函数，nm为节
点上的电子浓度；

[0064] 格边缘eij的SU电流密度模型表示为：

[0065]

[0066] 其中 是沿边缘eij的电场，μn,ij是边缘电子迁移率，Dn,ij是边缘扩散系数，是辅助人工扩散系数，其在算法性能方面起着关键作用；另外辅助人工扩散系数 满足
关系式

[0067]

[0068] 其中lij为边缘eij的长度，α为迎风系数， 为佩克莱特数，其中 为电势，μn以及Dn分别为载流子迁移率、载流子扩散率。

[0069] 优选的，所述电子电流密度计算模块具体如下：通过利用Nedelec边缘基空间将式(3)中的SU电流密度模型在 相邻段的中心进行插值，得到

[0070]

[0071] 其中 为电流密度， 为矢量基函数。

[0072] 优选的，单元矩阵方程的构造模块具体如下：

[0073] 通过将公式(2)、(5)代入公式(1)并在整个Ωe区域进行积分，得到：

[0074]

[0075] 其中h为纽曼边界条件；

[0076] 尽管公式(6)左边部分的第二项是在网格边缘所进行的操作，但是节点的未知数仍被计算，通过将边缘电流密度函数公式(3)代入公式(6)，因此公式(6)又写为：

[0077]

[0078] 其中当定点为i和j时，σr分别为‑1和1；

[0079] 通过遍历网格Ωe中所有顶点，得到公式(7)的矩阵形式：

[0080]

[0081] [Ke]是单位“迎风”刚度矩阵，其单元由公式(9)计算出

[0082]

[0083] [Me]是阻尼矩阵，其元素表示为：

[0084]

[0085] fe是列向量，其元素表示为：

[0086]

[0087] 在应用时域后向差分之后，公式(8)离散为：

[0088] ([Me]‑[Ke]){n}t+Δt＝[Me]{n}t+{fe}  (12)

[0089] 其中Δt是演化时间的步长。

[0090] 优选的，所述系统矩阵方程的构造模块具体如下：

[0091] 通过对求解域中的所有网格单元进行遍历，得出系统矩阵方程：

[0092] ([M]‑[K]){n}t+Δt＝[M]{n}t+{f}  (13)

[0093] 其中[K]是系统“迎风”刚度矩阵，[M]是系统阻尼矩阵；在稳态情况下，与时间有关的项均应该设置为0，得到：

[0094] ‑[K]{n}＝{f}  (14)

[0095] 公式(13)以及公式(14)两个矩阵方程分别用于求解瞬态情况和稳态情况。

[0096] 相较于FBSG方法，本发明对空间网格质量要求更低，从而使得收敛性更好；相较于SUPG方法，本发明在迎风函数的多维应用方面进行了优化；相较于FVFEM‑SG方法，本发明在
定义人工扩散系数方面更具灵活性。

附图说明

[0097] 图1是漂移扩散模型数值方案的流程图；

[0098] 图2是混合流线迎风有限元方法实施步骤流程图；

[0099] 图3是四面体网格单元示意图；

[0100] 图4是PN结在三个维度的结构以及掺杂密度分布图((a)一维(b)二维(c)三维)；

[0101] 图5是基于PN结结构的SU‑FEM与商业仿真软件结果对比图((a)一维(b)二维(c)三维)；

[0102] 图6是两种质量的二维PN结网格剖分图((a)高质量网格(b)低质量网格)；

[0103] 图7为SU‑FEM与SUPG仿真结果对比图(P端电压为0.2V时SU‑FEM与SUPG仿真结果对比图(a)电子密度(b)空穴密度)；

[0104] 图8为SU‑FEM与FBSG仿真结果对比图(P端电压为1.2V时SU‑FEM与FBSG仿真结果对比图(a)电子密度(b)空穴密度)；

[0105] 图9为SU‑FEM与CVFEM‑SG对半导体连续性方程离散程序示意图；

[0106] 图10为SU‑FEM与CVFEM‑SG对半导体连续性方程离散时间对比图；

[0107] 图11为混合流线迎风有限元系统框图。

具体实施方式

[0108] 以下通过特定的具体实施案例说明本发明的实施方式，本领域技术人员可由本说明书所描述内容直观地了解本发明的其他优点。本发明还可以通过另外不同的具体实施方
式加以实施或应用，本说明书中的各项细节也可以基于不同观点与应用，在没有背离本发
明的精神下进行各种修饰或改变。需说明的是，在不冲突的情况下，以下实施例及实施例中
的特征可以相互组合。

[0109] 本发明准确度验证：

[0110] 图1为漂移扩散模型数值方案的流程图，对漂移扩散模型的求解过程首先从生成网格、仿真设置以及初始化开始；其次在时间步范围内利用有限元方法对泊松方程进行数
值求解，并采用了牛顿‑拉弗森过程来提高数值计算的收敛性；然后，利用本发明所提出的
SU‑FEM对电子和空穴的电流连续性方程进行求解，以此得到电子和空穴的密度分布情况，
此过程具体实施方式对应图2中步骤2‑5，在下文详述。最后，当收敛条件满足时会开始新一
次的迭代过程。

[0111] 参照图2，本实施例用于半导体连续性方程的混合流线迎风有限元方法的具体步骤：

[0112] S1、几何模型的空间离散化：空间离散化可以将物理模型简化为每个网格上的时域差分方程，而在时域进行求解取决于器件的运行方式，在静态直流条件下，时间差可以直
接设置为零。

[0113] S2、网格单元的构造及应用：参照图3，其所展示的内容为顶点v1处的网格单元结构示意图，由图中可以看出每个网格单元由网格顶点、网格边缘、网格面组成。

[0114] 步骤2中，其关于电子连续性方程的积分过程具体如下：

[0115] S2.1.将电子连续性方程在网格单元上进行加权积分：

[0116]

[0117] 其中q为单位电荷，n为电子浓度，Jn为电子电流密度，Ωe为网格单元，Rn为净复合率，Nk为加权函数，dΩ为体积微元；

[0118] S2.2.通过标准拉格朗日插值函数可以得出每个网格内的电子浓度：

[0119]

[0120] m∈Ωe∪ΓN是指顶点vm位于计算域内或者纽曼边界上，m∈ΓD是指顶点vm属于狄利克雷边界，nd,m(t)表示相应的与时间相关的狄利克雷边界值，Nm为插值展开函数，nm为节
点上的电子浓度；

[0121] S2.3.为了更清楚地描述混合数值算法的构造，下面将只阐述在网格Ωe上的积分。为处理连续性方程的非线性，引入了SU电流密度模型来表征电子和空穴电流密度。沿网
格边缘eij的SU电流密度模型可以表示为：

[0122]

[0123] 其中 是沿边缘eij的电场，μn,ij是边缘电子迁移率，Dn,ij是边缘扩散系数，是辅助人工扩散系数，其在算法性能方面起着关键作用；另外辅助人工扩散系数 满足
关系式

[0124]

[0125] 其中lij为边缘eij的长度，α为迎风系数， 为佩克莱特数，为电势，μn以及Dn分别为载流子迁移率、载流子扩散率。

[0126] S3、电子电流密度计算：步骤3具体如下：

[0127] 混合流线迎风有限元方法通过利用Nedelec边缘基空间将式(3)中的SU电流密度模型在 相邻段的中心进行插值，因此可以得到

[0128]

[0129] 其中 为电流密度， 为矢量基函数；

[0130] S4、单元矩阵方程的构造：步骤4中将对单元矩阵方程的构造过程进行阐述：

[0131] S4.1.通过将公式(2)、(5)代入公式(1)并在整个Ωe区域进行积分，可以得到：

[0132]

[0133]

[0134] 其中h为纽曼边界条件；

[0135] S4.2.尽管公式(6)的左边部分的第二项是在网格边缘所进行的操作，但是节点的未知数仍被计算，通过将边缘电流密度函数公式(3)代入公式(6)，因此公式(6)又可写为：

[0136]

[0137] 其中当定点为i和j时，σr分别为‑1和1；

[0138] S4.3.通过遍历网格Ωe中所有顶点，可以得到公式(7)的矩阵形式：

[0139]

[0140] [Ke]是单位“迎风”刚度矩阵，其单元由公式(9)计算出

[0141]

[0142] [Me]是阻尼矩阵，其元素可以表示为：

[0143]

[0144] fe是列向量，其元素可以表示为：

[0145]

[0146] 在应用时域后向差分之后，公式(8)可以离散为：

[0147] ([Me]‑[Ke]){n}t+Δt＝[Me]{n}t+{fe}  (12)

[0148] 其中Δt是演化时间的步长；

[0149] S5、系统矩阵方程的构造：系统矩阵方程的形成具体步骤如下：

[0150] S5.1.通过对求解域中的所有网格单元进行遍历，可以得出系统矩阵方程：

[0151] ([M]‑[K]){n}t+Δt＝[M]{n}t+{f}  (13)

[0152] 其中[K]是系统“迎风”刚度矩阵，[M]是系统阻尼矩阵。另外，在稳态情况下，与时间有关的项均应该设置为0，因此可以得到：

[0153] ‑[K]{n}＝{f}  (14)

[0154] 公式(13)以及公式(14)两个矩阵方程可以分别用于求解瞬态情况和稳态情况。

[0155] 参照图4，图中展示了用来验证所提出数值方案准确度的三维PN结物理模型结构以及掺杂密度分布示意图，分别为每个掺杂区域定义了一个接触并且是具有狄利克雷边界
条件的纯欧姆型接触。从图5中可以看出，基于SU‑FEM的仿真器和COMSOL多物理场仿真软件
对PN结的仿真结果，二者具有很好的一致性，以此证明本发明所提出的数值计算方案在半
导体器件模拟中具有很好的准确性。

[0156] 上述实施例基于本发明所提出的用于漂移扩散模型的数值计算方案，并将该方案用于内部开发的仿真器，通过与COMSOL Multiphysics软件仿真结果对比后可以得出：本发
明可以很好的适用于半导体器件的仿真，其准确性已经达到了半导体器件仿真的要求。

[0157] 本发明性能验证：

[0158] 由上文所述，前文所提及的三种方法SU‑FEM、FBSG以及SUPG在一维情况下几乎没有区别，而在二维以及三维情况下，几种算法性能就会出现差别。在多维结构下，SU‑FEM与
传统的FBSG对于网格的要求会更低，并且相对于SUPG，其在多维应用中寻找最优迎风系数
会较容易。因此，为了验证本发明所提出方法的性能，基于图4(b)所示的二维PN结结构进行
仿真验证，且仿真过程需要求解的未知数为4141，得出结论如下：

[0159] 首先，在如图6(a)所示的高质量网格以及P端电压为0.2V时，对二维的PN结结构采用本发明所提出的SU‑FEM以及FBSG方法进行仿真，其结果如图7所示，从图中可以发现，在
高质量网格下，FBSG方法得出的结果误差仍然较大。

[0160] 另外，在基于图6(b)所示的低质量网格以及P端电压为1.2V时，对二维的PN结结构采用本发明所提出的SU‑FEM以及SUPG方法进行仿真，其结果如图8所示，从图中可以发现，
SUPG方法得出的结果出现了非物理震荡情况。

[0161] 最后，鉴于图9以及图10所示的SU‑FEM与CVFEM‑SG两种方法对半导体连续性方程离散程序过程对比以及时间对比图可以发现，由于CVFEM‑SG要形成有限体积单元，所以要
比SU‑FEM花费更多时间，从而展示了本发明所提出的方法的时间高效性。

[0162] 参见图11，一种用于半导体连续性方程的混合流线迎风有限元系统，其包括如下模块：

[0163] 几何模型的空间离散化模块：将物理模型简化为每个网格上的时域差分方程，而时域进行求解取决于器件的运行方式，在静态直流条件下，时间差设置为零。

[0164] 网格单元的构造及应用模块：所述网格单元为四面体形，由网格顶点、网格边缘以及网格面组成；网格单元的应用是将电子连续性方程在网格单元上进行积分；所述网格单
元的具体应用包括：

[0165] 将电子连续性方程在网格单元上进行加权积分：

[0166]

[0167] 其中q为单位电荷，n为电子浓度，Jn为电子电流密度，Ωe为网格单元，Rn为净复合率，Nk为加权函数，dΩ为体积微元； 为电子浓度随时间的变化率；

[0168] 通过标准拉格朗日插值函数得出每个网格内的电子密度：

[0169]

[0170] m∈Ωe∪ΓN是指顶点vm位于计算域内或者纽曼边界上，m∈ΓD是指顶点vm属于狄利克雷边界，nd,m(t)表示相应的与时间相关的狄利克雷边界值，Nm为插值展开函数，nm为节
点上的电子浓度；

[0171] 格边缘eij的SU电流密度模型表示为：

[0172]

[0173] 其中 是沿边缘eij的电场，μn,ij是边缘电子迁移率，Dn,ij是边缘扩散系数，是辅助人工扩散系数，其在算法性能方面起着关键作用；另外辅助人工扩散系数 满足
关系式

[0174]

[0175] 其中lij为边缘eij的长度，α为迎风系数， 为佩克莱特数，其中 为电势，μn以及Dn分别为载流子迁移率、载流子扩散率。

[0176] 电子电流密度计算模块：利用Nedelec边缘基空间将电流密度模型在网格单元相邻段的中心进行插值得到电流密度；所述电子电流密度计算模块具体如下：通过利用
Nedelec边缘基空间将式(3)中的SU电流密度模型在 相邻段的中心进行插值，得到

[0177]

[0178] 其中 为电流密度， 为矢量基函数。

[0179] 单元矩阵方程构造模块：构造单元矩阵方程；具体如下：

[0180] 通过将公式(2)、(5)代入公式(1)并在整个Ωe区域进行积分，得到：

[0181]

[0182] 其中h为纽曼边界条件；

[0183] 尽管公式(6)左边部分的第二项是在网格边缘所进行的操作，但是节点的未知数仍被计算，通过将边缘电流密度函数公式(3)代入公式(6)，因此公式(6)又写为：

[0184]

[0185] 其中当定点为i和j时，σr分别为‑1和1；

[0186] 通过遍历网格Ωe中所有顶点，得到公式(7)的矩阵形式：

[0187]

[0188] [Ke]是单位“迎风”刚度矩阵，其单元由公式(9)计算出

[0189]

[0190] [Me]是阻尼矩阵，其元素表示为：

[0191]

[0192] fe是列向量，其元素表示为：

[0193]

[0194] 在应用时域后向差分之后，公式(8)离散为：

[0195] ([Me]‑[Ke]){n}t+Δt＝[Me]{n}t+{fe}  (12)

[0196] 其中Δt是演化时间的步长。

[0197] 系统矩阵方程构造模块：通过对求解域中的所有网格单元进行遍历，得到系统矩阵方程，具体如下：

[0198] 通过对求解域中的所有网格单元进行遍历，得出系统矩阵方程：

[0199] ([M]‑[K]){n}t+Δt＝[M]{n}t+{f}  (13)

[0200] 其中[K]是系统“迎风”刚度矩阵，[M]是系统阻尼矩阵；在稳态情况下，与时间有关的项均应该设置为0，得到：

[0201] ‑[K]{n}＝{f}  (14)

[0202] 公式(13)以及公式(14)两个矩阵方程分别用于求解瞬态情况和稳态情况。

[0203] 本发明应用于半导体仿真领域中以对流为主导的非线性载流子输运方程的数值离散化方法——流线迎风有限元法(Streamline Up‑wind Finite Element Method,SU‑
FEM)，并具体描述了该混合数值方法的构造和分析过程。本发明能够较好地应用于半导体
器件领域，特别是针对漂移扩散模型中电流连续性方程的离散化。本发明所提出的SU‑FEM
方法相较于几种传统的数值离散方法有较明显的优势：相较于流线迎风Petrov Galerkin
方法，SU‑FEM在算法的多维应用方面进行了优化，使得在离散过程中寻找最优迎风系数更
加容易；相较于FBSG方法，对网格的质量要求更加宽松，收敛性更好；相较于有限体积有限
元(Finite Volume Finite Element Method,FVFEM)‑SG方法可以更加灵活的处理迎风函
数优化问题并且离散电流连续性方程需要更短的时间。因此，本发明方法可以较好的用于
半导体器件领域，并对其模拟仿真的发展起着推动作用。

[0204] 注意，上述仅为本发明的较佳实施例及所运用技术原理。本领域技术人员会理解，本发明不限于这里所述的特定实施例，对本领域技术人员来说能够进行各种明显的变化、
重新调整和替代而不会脱离本发明的保护范围。因此，虽然通过以上实施例对本发明进行
了较为详细的说明，但是本发明不仅仅限于以上实施例，在不脱离本发明构思的情况下，还
可以包括更多其他等效实施例，而本发明的范围由所附的权利要求范围决定。