Open Access

Numerical Simulation Analysis Mathematics of Fluid Mechanics for Semiconductor Circuit Breaker

 and    | Jul 27, 2021

Cite

Introduction

The drift-diffusion model is a physical model that describes the movement of carriers in semiconductor materials and has essential application backgrounds in the fields of remote sensing and solar power generation. The model consists of a Poisson equation and two carrier continuity equations. It is a system of tightly coupled and strongly nonlinear partial differential equations. Among them, the carrier continuity equation is a singular perturbation equation. Furthermore, it is understood that there is a sharp boundary layer and the value range is extremely large. This has brought great difficulties to its numerical solution [1]. Only by selecting appropriate variables and carefully designed discrete formats can good numerical simulation results be obtained. In addition, because the equations are strongly nonlinear, we need to develop targeted and efficient algorithms for solving nonlinear equations. Most of the open-source and commercial software used for semiconductor device simulation are serial or only support small-scale parallelism. Limited by its computational scale, it can only be used for numerical simulation of two-dimensional devices. There are few scalable algorithms in the literature that are suitable for large-scale parallelism. This paper focuses on the numerical calculation method of the three-dimensional drift-diffusion model. The purpose is to provide a scalable parallel algorithm suitable for large-scale equality for the numerical simulation of three-dimensional semiconductor devices.

The limited volume method is the discrete method used by most open-source and commercial semiconductor device simulation software. It has local conservation. The discrete system satisfies the principle of discrete extremum, which can keep the carrier concentration always positive. The discrete current flux in the finite volume method usually adopts the Scharfetter-Gummel format. It is an adaptive exponential upwind style, which adds artificial viscosity to the original problem and reduces numerical fluctuations. This paper uses this method as the basic numerical discrete format [2]. We use the coupled Newton iteration method to solve the nonlinear equations and propose using the algebraic multigrid (AMG) preconditioning GMRES method to solve the linear equations in the Newton iteration. Numerical experiments show that as long as multiple variables at the same geometric position are arranged together when sorting unknowns and treated as a whole when the grid is coarsened, the algebraic grid preconditioner can solve the problem efficiently. The number of iteration steps generally does not increase with the increase of the grid size, which ensures the parallel scalability of the algorithm. In the numerical experiment, the maximum grid size reached 500 million units, and the maximum number of processes reached 1024. This shows that our proposed algorithm has good scalability.

Drift-diffusion model

The steady-state drift-diffusion model consists of the Poisson equation −∇[(ɛψ) = q(pn + NAND), hole continuity equation ∇Jp = −qR, and electron continuity equation ∇Jn = −qR. The electrostatic potential ψ, the hole concentration p, and the electron concentration n are the physical quantities that need to be solved. ɛ is the dielectric constant of the material. q is the absolute value of the electron charge. NA and ND are the doping concentration of acceptor impurities and donor impurities, respectively. R is a compound term. Jp and Jn are hole current and electron current density, respectively. The Poisson equation is defined for the entire device (including semiconductors and oxides), and the continuity equation for electrons and holes is only defined for semiconductor materials. In the drift-diffusion model, Jp and Jn have the following forms, respectively Jp=qDpp+qμpp(ψ),Jn=qDnn+qμnn(ψ) \begin{equation}J_{p}=-q D_{p} \nabla p+q \mu_{p} p(-\nabla \psi), J_{n}=q D_{n} \nabla n+q \mu_{n} n(-\nabla \psi)\end{equation}

Here −qDpp and qDnn are the diffusion terms. The equation describes the diffusion motion of holes and electrons caused by the difference in concentration. p p(−∇ψ) and pn(−∇ψ) are drift terms [3]. The equation describes the drifting movement of holes and electrons generated by the electric field. Dp and Dn are the diffusion coefficients of holes and electrons, respectively. Dn is the hole mobility. μn is the electron mobility. Let k denote Boltzmann's constant. T represents temperature. Dp and μp, Dn, and μn satisfy the Einstein relationship: Dp=kTqμp,Dn=kTqμn \begin{equation}D_{p}=\frac{k T}{q} \mu_{p}, D_{n}=\frac{k T}{q} \mu_{n}\end{equation}

We introduce the quasi-Fermi levels φp and φn and the intrinsic carrier concentration nie, which are defined by the following relationship: p=nieexp(q(ψϕp)/kT,n=nieexp(q(ψϕn)/kT) \begin{equation}p=n_{i e} \exp \left(-q\left(\psi-\phi_{p}\right) / k T, n=n_{i e} \exp \left(q\left(\psi-\phi_{n}\right) / k T\right)\right.\end{equation}

Then the currents Jp and Jn can be written as Jp=qμppϕp,Jn=qμnnϕn \begin{equation}J_{p}=-q \mu_{p} p \nabla \phi_{p}, J_{n}=-q \mu_{n} n \nabla \phi_{n}\end{equation}

Common semiconductor device boundaries include ohmic contacts, metal-oxide limits, free borders, and so on.

Finite Volume Discrete
Selection of unknown variables

In addition to the electrostatic potential ψ, the unknown variables used in the numerical simulation of the drift-diffusion model can have different choices for the other two variables. This includes the carrier concentration (p, n) and the quasi-Fermi potential (φp, φn). In addition, we can also introduce Slotboom variables Φp and Φn: Φp=p/nieexp(qψ/kT),Φn=p/nieexp(qψ/kT) \begin{equation}\Phi_{p}=p / n_{i e} \exp (q \psi / k T), \Phi_{n}=p / n_{i e} \exp (-q \psi / k T)\end{equation}

In the actual calculation, one of the three groups of variables after the corresponding non-dimensionalization is often used: (ψ/kT,p/nie,n/nie),(ψ/kT,ϕp/kT,ϕn/kT),(ψ/kT,ϕp,ϕn) \begin{equation}\left(\psi / k T, p / n_{i e}, n / n_{i e}\right),\left(\psi / k T, \phi_{p} / k T, \phi_{n} / k T\right), \left(\psi / k T, \phi_{p}, \phi_{n}\right)\end{equation}

The variable group (ψ/kT, p/nie, n/nie) is the earliest used for numerical calculations. Because electric potential and carrier concentration are the physical quantities that people are most concerned about, it is most natural to choose this variable in physics. However, since the variation range of p and n reaches more than ten orders of magnitude, and the variation range of ψ/kT is only within two orders of magnitude, this set of variables is only suitable for decoupled Gummel iteration [4]. When using variable groups, ψ/kT, φp/kT and φn/kT have the same physical meaning and similar orders of magnitude. They are particularly ideal for coupled Newton iterations. The Jacobi matrix condition number obtained at this time is the smallest among the three sets of variables.

The variables Φp and Φn have no actual physical meaning. The range of changes of Φp and Φn is exceptionally drastic, so they are not suitable for use in precise calculations. However, they have essential value in theoretical analysis. Under this set of variables, the hole continuity equation and the electron continuity equation are transformed from asymmetric convective diffusion to a symmetric positive definite elliptic equation. On the one hand, this brings benefits to error analysis. On the other hand, some finite element methods for variable coefficient elliptic equations can be applied to semiconductor equations [5]. When constructing a numerical algorithm, you can first use this group of variables to build a discrete format and then replace it with (ψ/kT, p/nie, n/nie) or (ψ/kT, φp/kT, φn/kT) as the variable in the actual calculation.

Discrete finite volume

The carrier continuity equation in the drift-diffusion model is the convection-diffusion equation of strong convection. It is understood that there is a sharp boundary layer. The general finite element or finite difference method is easy to cause numerical oscillations. Therefore, we need to adopt the upside-down style or exponential fitting format. In addition, the ideal discrete structure should make the discrete system satisfy the principle of discrete extremum. It requires that the extreme value of the numerical solution is reached on the boundary and the carrier concentration is always positive [6]. This paper uses a three-dimensional tetrahedral unstructured grid. The control volume of the finite volume method is the area centered on the vertex and surrounded by the vertical bisector of each side. As shown in Figure 1, Figure 2.

Fig. 1

Control body centered on the vertex.

Fig. 2

Division of the control body in the unit.

Suppose xi is a vertex of the mesh, Ci represents the control body corresponding to the vertex, (xj| j = 1,2,, n) is the vertex adjacent to xi, and {lj| j = 1,2, , n} is the edge from xi. {Sj| j = 1,2, , n} is the boundary surface of the control volume bisected by lj. Note that Poisson equation (1) and carrier continuity equations (2), (3) can be written in the form of conservation laws as follows: F=R \begin{equation}\nabla F = R\end{equation}

Discrete it in a finite volume, and get the equation for each control body Ci: j=1nFj|Sj|=Ri|Ci| \begin{equation}\sum_{j = 1}^{n}F_{j}|S_{j}| = R_{i}|C_{i}|\end{equation}

Here, Fj represents the flux along the edge lj. The Poisson equation, Fj is given by Fj = ɛ(ψiψj)/|lj|. For the carrier continuity equation.

Scharfetter-Gummel format

The Scharfetter-Gummel format is a one-dimensional discrete format for calculating electron current and hole current density [7]. Define the Bernoulli function: B(x)=x/exp(x)1 \begin{equation}B(x) = x/\lbrack\exp(x) - 1\rbrack\end{equation}

The formula for calculating the electron current and gap current on the side in the SG format is jn|lj:=kTμnnie|lj|B(q(ψiψj)/kT)exp(qψi/kT)(ΦnjΦni)=kTμnnie|lj|B(q(ψiψj)/kT)exp(qψi/kT)(exp(qϕni/kT)exp(qϕnj/kT))=kTμn|lj|(njB(q(ψiψj)/kT)niB(q(ψiψj)/kT)) \begin{equation}\begin{aligned}&j_{n} \mid l_{j}:=\frac{k T \mu_{n} n_{i e}}{\left|l_{j}\right|} B\left(q\left(\psi_{i}-\psi_{j}\right) / k T\right) \exp \left(q \psi_{i} / k T\right)\left(\Phi_{n j}-\Phi_{n i}\right)= \\&\frac{k T \mu_{n} n_{i e}}{\left|l_{j}\right|} B\left(q\left(\psi_{i}-\psi_{j}\right) / k T\right) \exp \left(q \psi_{i} / k T\right)\left(\exp \left(-q \phi_{n i} / k T\right)-\exp \left(-q \phi_{n j} / k T\right)\right)= \\&\frac{k T \mu_{n}}{\left|l_{j}\right|}\left(n_{j} B\left(q\left(\psi_{i}-\psi_{j}\right) / k T\right)-n_{i} B\left(q\left(\psi_{i}-\psi_{j}\right) / k T\right)\right)\end{aligned}\end{equation} jp|lj:=kTμpnie|lj|B(q(ψjψi)/kT)exp(qψi/kT)(ΦpiΦpj)=kTμpnie|lj|B(q(ψjψi)/kT)exp(qψi/kT)(exp(qϕpi/kT)exp(qϕpj/kT))=kTμp|lj|(piB(q(ψjψi)/kT)pjB(q(ψiψj)/kT)) \begin{equation}\begin{aligned}&j_{p} \mid l_{j}:=\frac{k T \mu_{p} n_{i e}}{\left|l_{j}\right|} B\left(q\left(\psi_{j}-\psi_{i}\right) / k T\right) \exp \left(-q \psi_{i} / k T\right)\left(\Phi_{p i}-\Phi_{p j}\right)= \\&\frac{k T \mu_{p} n_{i e}}{\left|l_{j}\right|} B\left(q\left(\psi_{j}-\psi_{i}\right) / k T\right) \exp \left(-q \psi_{i} / k T\right)\left(\exp \left(-q \phi_{p i} / k T\right)-\exp \left(-q \phi_{p j} / k T\right)\right)= \\&\frac{k T \mu_{p}}{\left|l_{j}\right|}\left(p_{i} B\left(q\left(\psi_{j}-\psi_{i}\right) / k T\right)-p_{j} B\left(q\left(\psi_{i}-\psi_{j}\right) / k T\right)\right)\end{aligned}\end{equation}

The SG format is an adaptive exponential upwind style. When ψi = ψj, the electric field is zero, and there is only diffusion current inside the device without drift current [8]. At this time, the current given by the SG format is jn|lj=kTμn|lj|(njni),jp|lj=kTμp|lj|(pipj) \begin{equation}j_{n}|l_{j} = \frac{\text{kT}\mu_{n}}{|l_{j}|}(n_{j} - n_{i}),j_{p}|l_{j} = \frac{\text{kT}\mu_{p}}{|l_{j}|}(p_{i} - p_{j})\end{equation}

Corresponds to the central difference format. When ψiψj jn|lj=qμnψiψj|lj|ni,jp|lj=qμpψiψj|lj|pj \begin{equation}j_{n}|l_{j} = q\mu_{n}\frac{\psi_{i} - \psi_{j}}{|l_{j}|}n_{i},j_{p}|l_{j} = q\mu_{p}\frac{\psi_{i} - \psi_{j}}{|l_{j}|}p_{j}\end{equation} When ψiψj jn|lj=qμnψiψj|lj|nj,jp|lj=qμpψiψj|lj|pi \begin{equation}j_{n}|l_{j} = q\mu_{n}\frac{\psi_{i} - \psi_{j}}{|l_{j}|}n_{j},j_{p}|l_{j} = q\mu_{p}\frac{\psi_{i} - \psi_{j}}{|l_{j}|}p_{i}\end{equation}

In the above two cases, the current given by the SG format is equivalent to the first-order upwind formula.

Solving nonlinear equations
Nonlinear iterative method

The equation of the drift-diffusion model is {F1(ψ,ϕp,ϕn)=0F2(ψ,ϕp,ϕn)=0F2(ψ,ϕp,ϕn)=0 \begin{equation}\left\{ \begin{matrix}F_{1}(\psi,\phi_{p},\phi_{n}) = 0 \\F_{2}(\psi,\phi_{p},\phi_{n}) = 0 \\F_{2}(\psi,\phi_{p},\phi_{n}) = 0 \\\end{matrix} \right.\ \end{equation}

This is a system of nonlinear equations, which must be solved by an iterative method. Commonly used iterative methods include Gummel iteration for decoupled solution of equations and coupled Newton iteration.

Gummel iteration

The Gummel iteration method is a decoupled iteration method equivalent to the nonlinear block Gauss-Seidel iteration. It decouples the system of equations [9]. First fix φp and φn to solve Poisson equation (1), then select the updated electrostatic potential ψ, and solve equations (2) and (3) in turn. The iterative process is as follows: {F1(ψn+1,ϕpn,ϕnn)=0F2(ψn+1,ϕpn+1,ϕnn)=0F2(ψn+1,ϕpn+1,ϕnn+1)=0 \begin{equation}\left\{\begin{array}{l}F_{1}\left(\psi^{n+1}, \phi_{p}^{n}, \phi_{n}^{n}\right)=0 \\F_{2}\left(\psi^{n+1}, \phi_{p}^{n+1}, \phi_{n}^{n}\right)=0 \\F_{2}\left(\psi^{n+1}, \phi_{p}^{n+1}, \phi_{n}^{n+1}\right)=0\end{array}\right.\end{equation}

Newton iteration

The Newton iteration is a coupled iterative method. The Jacobi matrix of the equation system needs to be calculated in each iteration [10]. The iterative process is as follows: (F1/ψF1/ϕpF1/ϕnF2/ψF2/ϕpF2/ϕnF3/ψF3/ϕpF3/ϕn)(δψδϕpδϕn)=(F1nF2nF3n) \begin{equation}\left(\begin{array}{lll}\partial F_{1} / \partial \psi & \partial F_{1} / \partial \phi_{p} & \partial F_{1} / \partial \phi_{n} \\\partial F_{2} / \partial \psi & \partial F_{2} / \partial \phi_{p} & \partial F_{2} / \partial \phi_{n} \\ \partial F_{3} / \partial \psi & \partial F_{3} / \partial \phi_{p} & \partial F_{3} / \partial \phi_{n}\end{array}\right)\left(\begin{array}{l}\delta \psi \\ \delta \phi_{p} \\ \delta \phi_{n}\end{array}\right)=\left(\begin{array}{c}-F_{1}^{n} \\-F_{2}^{n} \\-F_{3}^{n}\end{array}\right)\end{equation} {ψn+1=ψn+δψϕpn+1=ϕpn+δϕpϕnn+1=ϕnn+δϕn \begin{equation}\left\{\begin{array}{l}\psi^{n+1}=\psi^{n}+\delta \psi \\ \phi_{p}^{n+1}=\phi_{p}^{n}+\delta \phi_{p} \\ \phi_{n}^{n+1}=\phi_{n}^{n}+\delta \phi_{n}\end{array}\right.\end{equation}

The advantage of Newton iteration is that it has a quadratic convergence rate, and the number of iteration steps is not sensitive to voltage. The disadvantage is that a larger-scale linear equation system needs to be solved. Some numerical methods are not easy to get an accurate Jacobi matrix. If an inexact Jacobi matrix is used, the convergence speed will be reduced to linear convergence [11]. Furthermore, Newton's iteration is sensitive to the selection of initial values. It is often calculated in sequence from low voltage to high voltage in a certain step size in device simulation. When solving, take the solution obtained from the last voltage as the initial value of the iteration. Therefore, the Newton iteration limits the step size of the voltage increase.

The solution of linear equations in Newton iteration

We use (ψ/kT, φp/kT, φn/kT) it as the method of solving the linear equations obtained by coupling Newton iteration with variables. The unraveling of linear equations in semiconductor numerical simulation is the most time-consuming part. The condition number of its coefficient matrix is vast, and the values of the matrix elements differ by more than ten orders of magnitude. For this kind of multi-physics model similar to the drift-diffusion model, the solution of linear equations obtained by coupled Newton iteration is still a hot research topic.

Most of the early methods used are sparse direct methods or block iterative approaches. The Jacobi matrix coupled with Newton iteration has a block structure. The block iteration method regards each matrix block as an element and performs Richardson, Jacobi, or Gauss-Seidel iteration on a large matrix. However, the article pointed out that various block iteration methods are not satisfactory for the semiconductor device after startup. The number of iterative convergence steps will increase as the voltage increases, and its behavior is similar to the decoupled Gummel iteration. After many experiments, we have found that the existing algebraic multigrid (AMG) solver is produced by solving the coupled Newton iteration Very effective preconditioner for linear equations. Using the GMRES method and combining algebraic multigrid preconditioners can achieve high computational efficiency and good parallel scalability [12]. We need to consecutively different number variables at the same geometric position. We deal with them as a whole in the mesh coarsening module of the Algebraic Multigrid Presolver. The following numerical experiments show that we use the GMRES method combined with the algebraic multigrid preconditioner to solve the linear equations generated by the coupled Newton iteration. The number of iteration steps does not generally increase with the increase of the grid size.

Implementation of the parallel program based on PNG platform

We implemented the algorithm introduced in this article based on the parallel adaptive meta-software platform PNG. PNG provides the management of parallel three-dimensional tetrahedral unstructured grids, supports large-scale parallelism, and has good scalability. Although it is designed to develop parallel adaptive finite element programs, it can also support the development of limited volume algorithm programs.

The Newton iterative algorithm in our program is directly implemented based on the PHG's degree of freedom and vector management function. The linear equations are solved by calling the GMRES solver in the parallel open-source solver PETSc to complete. The algebraic multigrid preconditioner uses Boomer AMG, another parallel open-source solver Hyper. We found in numerical experiments that such a combination has the best solution efficiency and robustness [13]. The coarsening method of Boomer AMG preconditioner is “Falgout,” which corresponds to the 3 × 3 sub-block structure of the linear equation system. We set the parameter of HYPRE_Boomer AMG Set Num Functions to 3.

Numerical examples

The calculation examples used in the thesis are PN junction and MOSFET transistors. The cluster contains 282 computing nodes, including two Intel X55504 core CPUs and 24GB of memory. The interconnection network is Infiniband DDR.

PN junction

The device structure and dimensions are shown in Figure 3. The doping concentration of the P area and N area is both 1 × 1017 cm−3. Figure 4 shows the IV curve of the PN junction. It can be seen from the numerical results that although the problem we are solving is a heavily doped semiconductor device. Carriers inside the device have powerful convective motion. Therefore, the calculation result has a very sharp boundary layer. However, the algorithm used in this paper has a good description of the boundary layer while avoiding numerical oscillations.

Fig. 3

PN junction.

Fig. 4

IV curve of PN junction.

MOSFET

This is an example of heavy doping concentration. The device structure and dimensions are shown in Figure 5. The thickness of the oxide layer is 0.025 nm, the doping concentration of the source region and the drain region are both 1 × 1020 cm−3 and the substrate doping concentration is 1 × 1016 cm−3.

Fig. 5

Metal-oxide-semiconductor field-effect transistor.

Figure 6 shows the output characteristic curve of the MOSFET. Compared with the PN junction, the carrier movement in the MOSFET is more complicated. The device can work in different modes such as the cut-off region, the saturation region, and the linear region. It can be seen from the calculation results that our algorithm is highly consistent with the theoretical results under these modes.

Fig. 6

MOSFET output characteristic curve.

Parallel scalability

The primary calculation time for semiconductor device simulation lies in the solution of linear equations. In large-scale parallel computing, the resolution of linear equations is based on iterative methods. Generally speaking, as the grid size increases, the number of iteration steps of the linear equation system increases. Table 1 and Table 2 respectively show the number of Newton iteration steps (NT) required for PN junction and MOSFET simulation under different grid sizes and the number of MPI processes, as well as the average number of iteration steps (GMRES) and average solution time (Time). It can be seen from the table that the number of iterative steps for solving the linear equations in the algorithm of this paper does not increase with the increase of the grid size and the number of processes, and it is not sensitive to the rise in voltage. This illustrates the efficiency of algebraic multigrid preconditioners.

Numerical simulation of PN junction diffusion test.

Voltage/V 8429568 unit 16 processes 67436544 unit 128 processes

NT GMRES Time/s NT GMRES Time/s
0.2 6 25 27.7 6 29 53.8
0.4 6 27.3 29.3 6 27.2 52
0.6 6 30.3 31.3 6 32 57.2
0.8 6 20.7 24.7 6 24 48.5
1 4 20.8 20.3 4 24.8 45.8
1.2 3 26 23.2 3 32 50.4
1.4 3 35.3 28.7 3 39.3 55.2
1.6 3 45.7 34.3 3 49.7 61.5
1.8 3 53.7 38.7 3 62.3 72.7
2 3 61.7 43.3 3 74 82.4

MOSFET numerical simulation scalability test (VG=8V).

VD/V 5050368 Unit 8 Process 40402944 unit 64 processes

NT GMRES Time/s NT GMRES Time/s
2 6 35.5 44.8 6 49.7 87.7
4 5 32.2 41.6 5 47.4 84.8
6 5 30.4 39.7 5 40 75.2
8 5 32 41.4 5 41 76.2
10 5 32 41.2 4 45.2 81.5
12 4 37 45.8 4 45.5 81.8
14 4 38.2 46.8 4 43 78.9
16 4 37.2 45.9 4 44.8 80.8
18 4 39 47.5 4 46 82.3
20 4 38.8 47.4 4 45.8 81.8
Conclusion

The paper proposes a parallel algorithm for the numerical simulation of the three-dimensional drift-diffusion model. This method uses the finite volume method for numerical discretization. In addition, it uses the coupled Newton iteration and the GMRES iteration based on algebraic multigrid preconditions to solve the nonlinear equations generated by the discretization. Numerical experiments show that the proposed method is robust, effective, and highly scalable and can be applied to three-dimensional large-scale parallel numerical simulations of some semiconductor devices.

eISSN:
2444-8656
Language:
English
Publication timeframe:
Volume Open
Journal Subjects:
Life Sciences, other, Mathematics, Applied Mathematics, General Mathematics, Physics