混合流线迎风有限体积有限元方法、模型数值离散系统

摘要

本发明公开了用于半导体连续性方程的混合流线迎风有限体积有限元方法、模型数值离散系统。本发明混合流线迎风有限体积有限元方法：几何模型的空间离散化；构造控制体积单元的；利用边缘电流密度和矢量基插值来计算出网格内部电流密度；构造单元矩阵方程的；构造系统矩阵方程的。本发明能够较好地应用于半导体器件中漂移扩散模型中电流连续性方程的离散化，对于半导体器件地模拟仿真发展有着推动作用，并且在固体电子领域的数值建模中，相较于流线迎风Petrov Galerkin方法，解决了多维问题求解中边缘计算精度差问题，收敛性更好；相较于FBSG方法，对空间网格质量要求更加宽松，收敛性更好；相较于传统控制体积有限元方法，在迎风函数的优化方面更具灵活性。

说明书

技术领域

本发明属于半导体仿真技术领域，涉及半导体器件中漂移扩散模型中电流连续性方程的离散化，具体是一种用于半导体连续性方程的混合流线迎风有限体积有限元方法及基于混合流线迎风有限体积有限元方法的离散漂移扩散模型数值系统。

背景技术

数值建模是一种实用的计算机辅助设计工具，用于模拟和数字应用的半导体器件的设计和优化。在固体电子领域，漂移扩散模型(drift-diffusion(DD)model)是最常用的数学物理模型。漂移扩散模型是基于电子和空穴电流连续性方程与泊松方程耦合的漂移扩散图。然而，电子和空穴电流连续性方程的混合性质，即扩散对流特征和强非线性，使它们的空间离散化至关重要。

最初，研究人员试图基于常规的有限差分方法(FDM)或有限元方法(FEM)方案来求解电流连续性方程，但是这些方法会产生非物理振荡。通常，可以通过对常规数值算法进行附加的人工扩散来实现稳定性。附加的辅助人工扩散率不仅需要稳定数值算法，而且还应避免混淆数值解。因此，需研发一种有效的人工扩散率的构造方法。目前广泛应用的FB-SG(finite boxes-Scharfetter-Gummel)算法对于网格要求严格，其需要通过控制体积有限元方案与边缘SG电流模型相结合来克服该缺点。流线型迎风Petrov Galerkin(SUPG)方法是离散电子和空穴电流连续性方程的另一种有效方法，然而此方法为多维应用程序寻找“最优”的迎风函数极具挑战性。

发明内容

为解决现有技术存在的上述问题，本发明基于控制体积有限元方案提出一种混合流线迎风有限体积有限元方法(SU-FVFEM)以及基于该方法的数值离散系统。根据本发明方法以及数值离散方法对半导体电学特性和电热耦合进行仿真验证可知，相较于FBSG方法，本发明对空间网格质量要求更加宽松，收敛性更好；相较于SUPG方法，本发明解决了多维问题求解中边缘计算精度差问题，收敛性更好；相较于传统控制体积有限元方法，本发明在迎风函数的优化方面更具灵活性。

本发明采用以下技术方案：

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

S1、几何模型的空间离散化；

S2、构造控制体积单元；

S3、利用边缘电流密度和矢量基插值计算出网格内部电流密度；

S4、构造单元矩阵方程：混合有限体积-有限元格式和时域后向差分；

S5、构造系统矩阵方程：遍历所有网格单元。

进一步地，步骤S1空间离散化的具体步骤：将物理模型简化为每个网格上的时域微分方程，时域解的要求取决于器件的运行方式，对于静态直流条件，可以将时差直接设置为零。

进一步地，步骤S2中，控制体积单元由属于八个不同网格的八个部分组成，每个部分包括网格顶点、网格重心、面重心和边缘重心。

进一步地，步骤S2中，控制体积单元的应用：

S2.1.将电子连续性方程集成在控制体积单元上：

其中，q为单位电荷，n为电子浓度，J

S2.2.通过在两边应用散度定理得到其“弱”形式：

其中，

S2.3.通过标准拉格朗日基可以得出每个网格内的电子密度：

其中，m∈Ω

S2.4.引入流线引入流线型迎风(SU)电流密度模型来处理电子和空穴电流连续性方程的非线性，沿网格边缘的SU电流密度模型可以表示为(对网格①(见图3)的积分进行阐述)：

其中，E

S2.5.将q

进一步地对步骤S3进行阐述：

S3.1.对于结构网格，边缘SU电流密度模型可以计算

其中，W

进一步地对步骤S4中单元矩阵方程的构造进行阐述：

S4.1.通过将公式(3)-(5)代入公式(4)得到①∩Ω

其中等效的人工扩散率用公式

S4.2.尽管公式(7)的左边部分的第二项是对边缘的操作，但是节点的未知数仍被计算，因此公式(7)左边的第二项可以简化为：

其中，l

S4.3.通过遍历网格①中所有顶点，可以将公式(8)写成矩阵形式：

K

M

f

在考虑时域后向差分之后，公式(9)可以离散为：

([M

其中，Δt是演化时间的步长。

进一步地对步骤S5中系统矩阵方程的构造进行阐述：

S5.1.通过遍历求解域中的所有网格单元，可以给出系统方程：

([M]-[K]){n}

其中，[K]是系统迎风刚度矩阵，[M]是系统阻尼矩阵，而对于稳态情况，时间依赖项设置为0得到：

-[K]{n}＝{f} (15)

式(13)和(15)两个矩阵可以分别用于求解瞬态情况和稳态情况。

本发明还公开了一种基于混合流线迎风有限体积有限元方法离散漂移扩散模型数值系统，其包括如下模块：

初始仿真设置模块：几何模型空间离散化，参数以及存储初始化；

泊松方程求解模块：使用有限元算法结合牛顿-拉弗森过程求解泊松方程得到电势分布；

电子和空穴电流连续性方程求解模块：将电势分布作为输入条件求解电子和空穴的电流连续性方程，使用上述的混合流线迎风有限体积有限元方法求解电子和空穴的电流连续性方程得到载流子浓度分布。

本发明的优点在于：

1.相较于FBSG方法，本发明对空间网格质量要求更加宽松，收敛性更好；

2.相较于SUPG方法，本发明解决了多维问题求解中边缘计算精度差问题，收敛性更好；

3.相较于传统控制体积有限元方法，本发明在迎风函数的优化方面更具灵活性。

附图说明

图1是基于混合流线迎风有限体积有限元方法的漂移扩散模型数值处理流程图；

图2是一种混合流线迎风有限体积有限元方法的离散过程图；

图3是在顶点ν

图4是具有不同迎风函数的数值算法的相对误差坐标图；

图5是三维PN结二极管的几何示意图和离子掺杂密度分布图；

图6是使用内部开发的模拟器和COMSOL多物理场模拟器对PN结二极管的模拟结果图；

图7是电势、电子密度和空穴密度的3-D分布图；

图8是电热耦合过程的流程图；

图9是基于混合流线迎风有限体积有限元方法的漂移扩散模型数值系统框图。

具体实施方式

以下结合附图详述本发明的优选实施例。本发明还可以通过另外不同的具体实施方式加以实施或应用，本说明书中的各项细节也可以基于不同观点与应用，在没有背离本发明的精神下进行各种修饰或改变。需说明的是，在不冲突的情况下，以下实施例及实施例中的特征可以相互组合。

实施例一：

图1是漂移扩散模型数值处理流程图，图9为系统框图，基于本发明提出的混合流线迎风有限体积有限元方法提出了一种数值方案，该方案使用有限元方法(FEM)离散泊松方程，用所提出的混合流线迎风有限体积有限元方法(SU-FVFEM)求解电子和空穴电流连续性方程，并引入牛顿-拉弗森(N-R)过程来提高数值计算的收敛性，从而实现漂移扩散模型数值求解，具体步骤如下：

初始仿真设置模块，用于初始仿真的设置：几何模型空间离散化，各项参数以及存储初始化；

泊松方程求解模块：使用有限元算法结合牛顿-拉弗森过程求解泊松方程得到电势分布；

电子和空穴电流连续性方程求解模块：将电势分布作为输入条件求解电子和空穴的电流连续性方程，使用SU-FVFEM求解电子和空穴的电流连续性方程得到载流子浓度分布(详见下文实施例叙述)；在此过程中判断是否有迭代溢出的现象，若无迭代溢出，判断是否收敛，若计算过程达到收敛条件，则输出结果；若未达到收敛条件，则重新开始泊松方程的求解；若有迭代溢出，则在输出结果中显示迭代溢出，计算中止。

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

S1、几何模型的空间离散化；

S2、控制体积单元的构造：包括网格顶点、网格重心、面重心和边缘重心；

S3、利用边缘电流密度和矢量基插值来计算出网格内部电流密度；

S4、单元矩阵方程的构造：混合有限体积-有限元格式和时域后向差分；

S5、系统矩阵方程的构造：遍历所有网格网格单元。

参照图3，顶点ν

步骤S2中，电子连续性方程在控制体积单元上的积分可以表现为以下步骤：

S2.1.将电子连续性方程集成在体积单元上：

其中，q为单位电荷，n为电子浓度，J

S2.2.通过在两边应用散度定理得到其“弱”形式：

其中，

S2.3.通过标准拉格朗日基可以得出每个网格内的电子密度：

其中，m∈Ω

S2.4.引入流线引入流线型迎风(SU)电流密度模型来处理电子和空穴电流连续性方程的非线性，沿网格边缘的SU电流密度模型可以表示为(对网格①的积分进行阐述)：

其中，E

步骤2.4中，α在算法性能中起决定性作用并且为SU电流密度模型选好一个“最优”的α是一项关键任务，下面表1列出了四种关于α的定义式：

表1迎风系数α的定义

利用给出的数值方案对PIN结半导体进行模拟仿真，找出四项里最优的α参数。参照图4给出了使用不同迎风函数情况的相对误差，而第二和第三种迎风函数的定义在收敛方面的性能优于其他定义，在下面的工作中选取第三种迎风函数的定义，因其定义相较于第二种较为简单；

S2.5.将q

步骤S3具体阐述如下：

S3.1.对于结构网格，边缘SU电流密度模型可以计算

其中，W

对步骤S4中单元矩阵方程的构造进行阐述：

S4.1.通过将公式(3)-(5)代入公式(4)得到①∩Ω

其中等效的人工扩散率用公式

S4.2.尽管公式(7)的左边部分的第二项是对边缘的操作，但是节点的未知数仍被计算，因此公式(7)左边的第二项可以简化为：

其中，l

S4.3.通过遍历网格①中所有顶点，可以将公式(8)写成矩阵形式：

K

M

f

在考虑时域后向差分之后，公式(9)可以离散为：

([M

其中，Δt是演化时间的步长。

对步骤S5中系统矩阵方程的形成进行阐述：

S5.1.通过遍历求解域中的所有网格单元，可以给出系统方程：

([M]-[K]){n}

其中，[K]是系统“迎风”刚度矩阵，[M]是系统阻尼矩阵，而对于稳态情况，时间依赖项设置为0得到：

-[K]{n}＝{f} (15)

式(13)和(15)两个矩阵可以分别用于求解瞬态情况和稳态情况。

参照图5，展示了用来验证所提出数值方案性能的3D PN结二极管几何示意图和离子掺杂密度分布，该PN结二极管的尺寸是0.4μm×0.4μm×0.8μm。

参照图6可以得到基于该数值算法的内部开发的仿真器和COMSOL多物理场仿真器对PN结二极管的模拟结果。在模拟仿真的过程中，底部电极接地，顶部电极的电压为0-1V。图6(a)显示了两者的I-V特性曲线具有良好的一致性，图6(b),6(c),6(d)中给出了沿垂直轴对称线电势，电子，空穴的分布情况。

参照图7，COMSOL Multiphysics软件(图7(a))和基于该数值算法开发的内部仿真器(图7(b))实现的电势，电子和空穴的3D分布也很好地匹配，即开发的数值方案被证明能够在半导体器件模拟。

本发明提出一种离散漂移扩散模型的数值方案，并将该方案用于内部开发的仿真器，通过与COMSOL Multiphysics软件对比仿真后可以得出：本发明可以很好的用于半导体器件的模拟，并且在固体电子领域的数值建模中，相较于SUPG方法，本发明解决了多维问题求解中边缘计算精度差问题，收敛性更好；相较于FBSG方法，本发明对空间网格质量要求更加宽松，收敛性更好；相较于传统控制体积有限元方法，本发明在迎风函数的优化方面更具灵活性。

实施例二：

参照图8，本发明用于半导体的电热耦合仿真。在该数值方案中，分别利用SU-FVFEM，N-R(FEM)和时域有限元(FETD)方法离散电流连续性方程，泊松方程和热传导方程，电热耦合过程的耦合路径包括直接耦合(自热效应)和间接耦合(非线性材料参数)。主要步骤如下：

使用SU-FVFEM和N-R(FEM)求解漂移扩散模型，输出电势分布和载流子浓度分布；

基于得到的电势和载流子浓度分布，求得器件内部功率密度分布，更新非线性材料参数；

以求得的功率为输入求解热传导方程，输出温度分布；

基于电场分布和温度分布，更新非线性材料参数，并将其作为输入条件求解漂移扩散模型，直至耦合过程收敛。

本发明能够较好地应用于半导体器件中漂移扩散模型中电流连续性方程的离散化，对于半导体器件地模拟仿真发展有着推动作用，并且本方法相较于FBSG方法，对空间网格质量要求更加宽松，收敛性更好；相较于SUPG方法，解决了多维问题求解中边缘计算精度差问题，收敛性更好；相较于传统控制体积有限元方法在迎风函数的优化方面更具灵活性。

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