Numerical Simulation Analysis Mathematics of Fluid Mechanics for Semiconductor Circuit Breaker

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: (5)Φ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: (6)(ψ/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/n_ie, n/n_ie) 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/n_ie, n/n_ie) or (ψ/kT, φ_p/kT, φ_n/kT) as the variable in the actual calculation.

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

Fig. 2

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

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

Here, F_j represents the flux along the edge l_j. The Poisson equation, F_j is given by F_j = ɛ(ψ_i − ψ_j)/|l_j|. For the carrier continuity equation.

The Scharfetter-Gummel format is a one-dimensional discrete format for calculating electron current and hole current density [7]. Define the Bernoulli function: (9)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 (10)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}(11)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 (12)jn|lj=kTμn|lj|(nj−ni),jp|lj=kTμp|lj|(pi−pj)\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(13)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(14)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.

The equation of the drift-diffusion model is (15){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.

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.

The device structure and dimensions are shown in Figure 3. The doping concentration of the P area and N area is both 1 × 10¹⁷cm⁻³. 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

Fig. 4

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 × 10²⁰cm⁻³ and the substrate doping concentration is 1 × 10¹⁶cm⁻³.

Fig. 5

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

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.