一种描述硅PIN二极管反向恢复过程的二阶有限元方法

翟 弋, 张满红， 宫婷婷

(华北电力大学 电气与电子工程学院，北京 102206)

0 引 言

PIN二极管作为一种重要的功率器件，既可以单独被使用，也可以作为IGBT器件的一个元件被使用[1-2].它们通常在高注入水平下工作，其中具有运动边界的双极扩散方程(ADE，ambipolar diffusion equation)可以用来模拟器件的稳态和瞬态行为[3-4].现在研究人员已经提出了几种求解ADE的方法，如傅里叶展开法(FE)[5-6]、有限差分法(FD)[7-9]、有限元法(FEM)[10-11]和数值法[12].有限差分法和有限元法通常在SPICE模拟器中实现，而傅里叶展开法通常在MATLAB Simulink中实现[5].上述方法中，每种方法都有不同的近似值.

在已报道的文献[9]中，已经比较了描述硅PIN二极管反向恢复过程的有限元法和有限差分法.研究发现，在某些情况下，当空间格点的位置和相邻距离在迭代过程中被更新，通过边界载流子密度的零值位置来确定未耗尽N-区域2个运动边界的新坐标时，应用传统的FD方法所得到的p(x,t)的数据是完全错误的[7-8].为了解决这一问题，本文中，笔者提出了一种改进的FD方法，即利用单步后退欧拉法求解ADE时，在一次迭代中固定空间离散点.在每个时间步长收敛后，基于新的边界坐标建立新的空间网格，并用3次样条插值将p(x,t)从旧的空间网格转移到新的空间网格.按照这种方法，从FE法和改进的FD法中得到的具有相同参数集的p(x,t)数据与Silvaco Atlas模拟结果吻合良好[9].

通过Matlab仿真，笔者比较了用FE法和FEM法描述硅PIN二极管反向恢复过程的精度.在文献[10-11]中，将求解ADE问题转化为变分方程的极小化问题，可用FEM法求解.研究表明，为从变分函数中推导出ADE，空间求解区域的大小应是固定的，并且应用FEM基函数尽可能精确地估计∂p/∂t的值.为说明在FEM方法中改变空间求解区域的影响，将FEM法和FE基函数相结合，得到了一组新的载流子浓度方程，与原来的FE方程相比少了一项，实现了单一有限元.通过比较2组方程的解，检验不同空间求解区域的影响.

对于FEM方法中的多重有限元实现，为使p(x,t)的FEM基函数展开连续且可微，笔者提出了一种二阶有限元方法(SFEM)，并在描述硅PIN二极管反向恢复过程时将其性能与线性有限元法(LFEM)进行了比较.

1 理论模型

1.1 双极扩散方程

图1为PIN二极管和相应掺杂分布的示意图.N-区的宽度为W，从P+N-结(x=0)延伸到N-N+(x=W).在高注入条件下，载流子分布p(x,t)受以下ADE[9]控制：

图1 反向恢复过程中PIN二极管掺杂水平和载流子分布示意图

(1)

其中，Da=2DnDp/(Dn+Dp)表示双极性扩散系数，τhl表示大注入载流子寿命.xl表示P+N-结(表示为J1)耗尽区相对于x=0的厚度，并且正方向指向右侧.xr表示N-N+结耗尽区的厚度(表示为J2).(1)定义在未耗尽N-基区，范围从x=xl到x=W-xr.

在x=xl和x=W-xr的边界条件下，(1)可以写成[9-10]

(2)

在(2)中，In为电子电流，Ip为空穴电流，Idisp为J1和J2结耗尽区的位移电流.A为器件面积，q为电子电荷.Dn和Dp分别为电子和空穴的扩散系数.当xl>0时，忽略漂移扩散电子流，当xr>0时，忽略漂移扩散空穴电流.在ADE方法中无法确定这些参数.由于系统误差和数值误差，最终的载流子浓度分布p(x,t)并不完全符合(2).因此，在(3)中，xl和x=W-xr的2个新的空间偏导数也是根据三点FD近似计算出来的.在以下的讨论中，将(2)中的∂p/∂x表示为∂p/∂x|I，将(3)中的∂p/∂x表示为∂p/∂x|FD.

(3)

PIN二极管两端的总电压为

VD=IDRD+VJ1+VJ2-Vd1-Vd2.

(4)

第1项是未耗尽N-区的欧姆压降.VJ1和-Vd1分别是在正向偏压和反向偏压下J1结两端的电压.同样，VJ2和-Vd2是J2结的相对应电压值.对于Vd1，可以采用常用的Vd1∝(xl)2近似，其中所有可运动空穴在耗尽区以饱和速度移动，或者采用包含全空穴漂移速度—电场关系的模型.将第1种称为x2相关模型，将第2种称为全场模型.当未耗尽N-区宽度较小且2个结耗尽区接近穿通时，第2种模型效果更好.对于Vd2，很难解析推导出包含全场效应的表达式，因此只采用x2相关模型Vd2∝(xr)2，详见文献[9].

Silvaco Atlas TCAD和MATLAB仿真基于图2所示的电路.其中，I1为电流源，I1=50 A，Vs为电压源，Vs=200 V，L1为电感，R1为固定电阻，R1=1 mΩ，R2为可变电阻，其初始阻值为1 MΩ，以20 ns的时间常数按照指数衰减降低到1 mΩ.总电流ID的微分方程为

图2 Silvaco Atlas TCAD器件混合仿真电路图

I1R2(t)-(VD1+Vs)-(R1+RD(t)+R2(t))ID(t).

(5)

为了评估不同ADE求解方法中p(x,t)、ID及其分量的精度，用2种方法比较了N-区空间积分表面的载流子浓度.第1种方法中，在tλ时，N-区的空间积分表面载流子浓度可以表示为

(6)

为了提高上述积分的数值精度，利用3次样条插值方法，使用均匀离散点p(xi，tλ)上的数值来计算函数p(x，tλ)，同时，根据载流子数守恒的连续性方程，可以得到以下方程：

(7)

其中Idisp表示J2结的位移电流.原则上，(6)和(7)中Nt应该相同，但是，由于各种算法都会存在误差，它们可能会有所不同.为了方便，将(6)中的Nt表示为Npt，将(7)中的Nt表示为NIt，2个量之间的差异反映了载流子浓度是如何保持守恒的.

所有方程均采用隐式欧拉(BE)方法求解，则(7)可以改写为

(8)

1.2 有限元法

Buiatti等[8-9]使用以下函数建立FEM方法:

(9)

其中，第2行的表达式只针对一维情况.固定边界xb和xe的Jn和Jp，让p变为p+δp，则可以得到δp的一阶函数的变化，如下所示：

(10)

如果(p/D)∂(δp)/∂t这一项可以取0，则可以得到ADE和相应的边界条件.但是，上述变分函数不是关于时间的积分.在(9)中，∂p/∂t是前缀函数，不参与(10)中的变化，所以在(10)中不会出现(p/D)∂(δp)∂t这一项.因此，p(x,t)和∂p/∂t的确定就形成了一个自洽问题.FEM法的精度还取决于如何从FEM基函数近似计算∂p/∂t.此外，上述FEM方法是在假设空间求解区域大小固定的情况下推导出来的，即xb和xe的值不变.

在线性有限元法、线性一维基函数和M有限元中，对函数Π最小化求解，得到如下常微分方程组：

(11)

其中[M],[G],[L]均为矩阵形式，

(12)

(13)

(14)

Δx是每个元素的宽度.M和[M]是2个不同的量,前者表示未耗尽N-区中空间点的个数，后者表示(12)中的矩阵.

1.3 基于FE方法的有限元FEM方法

利用类似SPICE的模拟仿真工具，Chibante等[10-11]已经实现了上述的FEM方法,模拟仿真过程中，将硅PIN二极管中未耗尽N-区离散成10多个块，并将反向恢复过程的动态变化转化为小尺寸变化块的时间演化.以下的讨论中，使用FE方法来说明空间求解区域的大小变化对最终FEM结果精度的影响.

在FE方法中，用xb≡xl和xe≡W-xr表示N-区准中性部分随时间变化的左边界(开始)和右边界(结束).用空间坐标x-xb的傅里叶级数表示随时间变化的载流子浓度p(x,t)，其系数pk(t)随时间坐标t变化

(15)

(16)

在(15)中，傅里叶级数被截断，只保留M项.同时，Da和τhl是常数，在这里不考虑它们的空间相关性.边界条件是载流子浓度梯度∂p/∂x|xb和∂p/∂x|xe，以及2个边界速度dx/dt|xb和dx/dt|xe.

已知pk(t)，实空间点x1=xb，x2=xb+Δx，…，xM=xe的p(x,t)值可由以下矩阵向量积求得，其中Δx=(xe-xb)/(M-1)，

(17)

傅里叶展开式系数pk(t)应满足以下微分方程：

(18)

文献[9]已经求解了这组载流子浓度方程，这2个边界由反馈进行控制.将FE方法与FEM方法相结合，研究了在FEM方法中空间求解区域大小变化的影响.将整个未耗尽N-区看作一个单一有限元，并假设

(19)

(20)

用(15)表示p(x,t)，(19)表示∂p/∂t，(20)表示δp(x,t)以及以下假设：

(21)

函数Π的最小化得到pk(t)的另外一个方程组.新的方程组与方程组(18)几乎相同，但(18)的第2式中没有以下项：

(22)

通过(21)，推导过程中仍然假设∂p/∂t是准确的.因此，在有限元FEM实现中，(22)中的缺失项是由运动边界xb和xe引起的空间求解区域的大小变化而引起的.如果边界变化不是很快，那么最终的FEM载流子浓度方程仍然有很好的近似，是可用的.因此可以使用FE方法模拟分别包含和不包含(22)这一项的PIN二极管的反向恢复过程，以检验在FEM方法实现中尺寸变化的空间求解区域对最终结果的影响.将这种方法称为FE-FEM方法，一种基于FE方法的FEM方法，这种FEM方法只有一个单一有限元.

1.4 二阶有限元法

LFEM方法中，线性基函数的展开使得p(x,t)和∂p/∂t空间连续，但在空间节点处不可微，这将降低LFEM方法的精度.因此探讨一种二阶有限元方法(SFEM),通过比较Silvaco Atlas模拟仿真结果和FEM方法的计算结果，来评估FEM方法的精度.

在SFEM方法中，对于归一化坐标ξ∈[0，1]有4个基函数，如下所示：

(23)

可以将这4个形函数及其导数改写为

(24)

定义如下2个4×4矩阵[A1]和[B1]：

(25)

对于单个元素，x∈[x1，x2]，x2=x1+Δx，用4个形函数展开p(x,t)

(26)

在x1和x2处，p(x,t)对x的导数分别为α1/Δx和α2/Δx.函数Π对αi的最小化得到以下的有限元FEM方程:

(27)

将上述方程扩展为由M个等距点[x1，x2，…，xM]组成的多个有限元，得到载流子浓度p(x,t)的展开式，如下所示：

(28)

其中所有形函数φi的参数都限制在[0,1].对于i=1，2，…，2M，αi的最终FEM方程为

(29)

从矩阵[A1]和[B1]可以导出矩阵[A2]和[B2].将(29)写成2M×2M的矩阵形式

(30)

(31)

其中，α=[α1，α2，…，α2M]T，α1，α3，…，α2M-1给出了p(x,t)在空间点x1，x2，…，xM的值，α2/Δx，α4/Δx，…，α2M/Δx给出了∂p(x,t)/∂x在上述空间点的值.

Δx(tλ+1)=[W-xr(tλ+1)|s-xl(tλ+1)|s]/(M-1).

(32)

根据插值p(x,t)可以得到新的空间点集所对应的p(x,t)值，然后开始下一个时间步长的迭代.

除了3次样条插值外，线性插值也是可能的.在时间tλ+1时刻，可以得到p(xi|λ，tλ+1)，并且能够通过(33)得到近似的p(xi|λ+1，tλ+1).

(33)

(33)的第1个方程式给出了新坐标中第i个空间点相对于旧坐标的偏移量.在实际计算中，用适当的有限差分近似来代替(33)第2个方程式中的偏导数.与3次样条插值相比，线性插值的计算速度更快.试验表明，对于较大的M，这2种插值方法产生的结果几乎相同；对于较小的M值，如M=4，3次样条插值的结果稍好一些.

1.5 仿真参数

正向偏压稳态下的初始载流子浓度为

(34)

对于Silvaco Atlas模拟仿真[14]，包括了默认的载流子浓度相关迁移率模型conmob、默认的平行电场相关迁移率模型fldmob和默认的载流子浓度相关的SRH复合模型consrh以及默认的高掺杂诱导带隙窄化模型bgn.整个结构采用了2组consrh寿命参数τn0和τp0(慢和快).τn和τp是3个区中每个区的实际SRH寿命参数.在N-层中，ADE(1)中的参数τhl=τn+τp.考虑带隙收缩效应，根据文献[3-4]的方程计算理论hp和hn参数.所有的模拟仿真参数如表1所示.We和Nae分别是P+接触层的宽度及其掺杂水平.Ws和Nds是N+接触层的宽度及其掺杂水平.ΔEge是P+接触层相对于N-层的带隙收缩量，ΔEgs是N+接触层相对于N-层的带隙收缩量.

表1 模型仿真参数值

利用MATLAB计算基于ADE的反向恢复过程.为避免在数值计算中出现下溢或上溢现象，引入以下基本标度：长度l0=10-4cm=1 μm，时间t0=10-12s=1 ps，电压V0=1 V，电荷q0= 10-12C=1 pC，其他参量可以很容易地推导出.在计算过程中，时间步长从0.01 ps改为1 ps.

在之前的研究中已经对采用FE，FD和Silvaco Atlas 3种方法描述硅PIN二极管的反向恢复过程作了比较.在慢变情况下，FE和FD法在x2相关和全场Vd1模型中的载流子浓度分布p(x,t)仅与Silvaco Atlas的载流子浓度分布近似一致，这是因为P+N-结耗尽区的电子浓度很高，足以影响Silvaco Atlas模拟的耗尽宽度(xl)，但所有基于ADE的方法都忽略了这一点，从而给出了稍错误的xl和xr值.在快变的情况下，前期FE法和FD法中的x2相关和全场Vd1模型的p(x,t)曲线与Silvaco Atlas模拟的p(x,t)曲线吻合较好；而在后期，只有FE法和FD法的全场Vd1模型的p(x,t)曲线与Silvaco Atlas模拟的p(x,t)曲线吻合较好.x2相关的Vd1模型导致在J1结和J2结2个耗尽区之间出现穿通现象.仍然使用2组参数来检验FEM方法在描述硅PIN二极管反向恢复过程中的性能，只考虑全场Vd1模型.FD方法所得结果与FE方法所得结果相似，不作讨论.

2 结果与讨论

图3a给出了采用全FE法和FE-FEM法的慢变和快变的I/V曲线，其中不包含式(22).未耗尽N-区的空间点总数为25.图3b给出了全FE法和FE-FEM法对2个结耗尽宽度xl和xr的影响.可观察到xr明显增加，但幅度不大.结果表明，在FEM方法中，采用一个大小缓变的空间求解区域是一种很好的近似方法.

solid.FE; dot.FE-FEM; a.M=25时慢变和快变的反向恢复I/V曲线；b.J1结和J2结耗尽宽度xl和xr.

图4和图5中以FE法作为参考来检验LFEM法和SFEM法的性能.首先比较了M=25时FE法、LFEM法和SFEM法的结果.图4a,b分别给出了在FE法和LFEM法中，慢变和快变2种情况分别在x=xl，W-xr边界处Δ(dp/dx)=∂p/∂x|I-∂p/∂x|FD的差异.SFEM方法中，在x=xl，W-xr边界分别使用α2/Δx和α2M/Δx代替∂p/∂x|FD，在SFEM方法中，Δ(dp/dx)的幅值最小.FE法中，未耗尽N-区的边界是通过经验反馈控制的，p(xl,t)和p(W-xr，t)的值不大，通常也不接近于0.因此,∂p/∂x|FD只是恰好遵循了∂p/∂x|I的形状，但这两者有很大的区别.相比之下，由于采用了p(x,t)的连续可微的FEM基函数展开，SFEM方法中Δ(dp/dx)很小.LFEM法的性能介于FE法和SFEM法之间.FE方法中边界载流子浓度的松散控制也会影响载流子守恒.图5a，b分别给出了3种方法在M=25时ΔN=Npt-NIt的差异，ΔN在SFEM方法中同样是最小的.FE法和FEM法的性能差异是由于未耗尽N-区边界的反馈控制不能很好地保持载流子守恒造成的，而SFEM在这方面最好.

图4 FE法、LFEM法和SFEM法中Δ(d p/d x)=∂ p/∂ x|I-∂ p/∂ x|FD的比较

图5 FE法、LFEM法和SFEM法中ΔN=Npt-NIt的比较

通过改变空间点的数目来比较LFEM和SFEM方法的性能.图6a,b分别比较了LFEM和SFEM方法中M=4，8，25时Δ(dp/dx)的差异，SFEM的性能比LFEM好.图7中M=4，8，25时ΔN差异也有类似的趋势.可以看出，对于LFEM和SFEM，M=4的数据结果比M=8，25的数据结果差.

图6 M=4,8,25时,快变的Δ(d p/d x)=∂ p/∂ x|I-∂ p/∂ x|FD的比较

图7 M=4，8，25时,快变的ΔN=Npt-NIt的对比结果

图8a,b，c分别给出了对于快变，M=4，8，25时，Silvaco Atlas(直线)、LFEM(空心符号)和SFEM(实心符号)的多个时间点的载流子分布p(x,t).由图8a，在LFEM方法中观察到p(x,t)有明显的峰值过冲.当M=4时，LFEM和SFEM方法中的xl与Silvaco Atlas模拟仿真数据不一致.图8b中，LFEM和SFEM方法中M=8时对应的数据具有更好的p(x,t)曲线，但仍与Silvaco Atlas模拟仿真数据存在差异(曲线c，d，e).图8c中，利用LFEM和SFEM方法得到的p(x,t)的数据在整个演化过程中与Silvaco Atlas模拟仿真的数据非常吻合，并且图8c表明SFEM方法与LFEM方法相比性能更好.

a.initial; b.0.274 μs; c.0.282 μs; d.0.286 μs; e.0.297 μs; a.M=4; b.M=8; c.M=25.

从图6～8的数据中可以得出结论，基于FEM基函数的较好的∂p/∂t近似，如采用SFEM方法或采用更多的空间点，可提高最终结果的精度.

LFEM方法和SFEM方法有另一个不同之处.SFEM方法中p(x,t)的基函数展开可以用于一次时间步长迭代后外推新空间点网格上p(x,t)的值，而LFEM方法却不能做到，这是因为后者只是连续而不可微的.由于使用了基函数式(23)，目前SFEM方法用于均匀的空间点网格，并可以扩展到非均匀的空间点网格.

据文献报道，集总电荷模型[15]和LFEM模型[10-11]已被成功地应用于SPICE模拟工具，对具有不同空间求解区域的IGBT器件进行建模.已经测试了现有的二极管偏置参数，并且通过改变空间步长来求解FEM方程[7-8],结果与Silvaco Atlas模拟结果不符[9]，据推测原因是由于模拟过程太快所致，如电阻R2的时间常数只有20 ns，后续应对此做深入研究.

3 结 论

综上所述，对有限元FEM方法的研究表明，为了能够从变分函数中推导出ADE，空间求解区域的大小应该是固定的，并且∂p(x,t)/∂t应使用FEM基函数尽可能精确地逼近并自洽确定.将FE方法与FEM方法相结合，形成一种新的FEM方法(FE-FEM)，得到了一组新的载流子浓度方程，其与原来的FEM方法相比少了一项，这一项与未耗尽N-区的2个边界的移动有关，因此在空间大小变化的求解区域中采用FEM方法可能会引入误差.通过比较FE法和FE-FEM法的计算结果分析了缺失项的影响,并在此基础上提出了一种二阶有限元法(SFEM)，该方法中p(x,t)的FEM基函数展开既连续又可微，将其与Silvaco Atlas模拟仿真、FE法和线性有限元法(LFEM)进行比较发现，SFEM法优于LFEM法.