3.1 Die Diskretisierung der Wärmeflußgleichung  

Der erste Teil des Abschnittes beschreibt die Diskretisierung der stationären Wärmeflußgleichung (2.68). Im zweiten Teil wird auf die implementierte Berechnung der transienten Wärmeflußgleichung eingegangen.

Die Implementierung der Wärmeflußgleichung in MINIMOS-NT erfolgt durch die Berechnung der Wärmestromdichte $\vec{S}_L$. Ausgangspunkt der eindimensionalen Diskretisierung ist die eindimensionale stationäre Gleichung für den Wärmestrom (2.64)

$$\begin{eqnarray} \vec{S}_L(T_L)=-\kappa_L(T_L) \cdot \frac{\mathrm{d}T_L}{\mathrm{d} x}\; . \end{eqnarray}$$ (3.1)

Die temperaturabhängige Wärmeleitfähigkeit $\kappa(T_L)$ kann mit unterschiedlichen Ansätzen modelliert werden. Die gewählten Funktionen sollten aufgrund der Differentialgleichung erster Ordnung (3.1) jedoch einmal analytisch integrierbar sein. In einer sehr allgemeinen Formulierung kann man $\kappa(T_L)$ durch ein Polynom ausdrücken

$$\begin{eqnarray} \kappa(T_L)=C_1+C_2\cdot (T_L-T_0)+C_3 \cdot(T_L-T_0)^2\; . \end{eqnarray}$$ (3.2)

Dies entspricht einer numerischen Taylorreihe um den Entwicklungspunkt $T_0\!\!=\!\!300\mathrm{K}$. Die Leitfähigkeit ist für die meisten Stoffe im entsprechenden Temperaturbereich gut bekannt. Daher ist es leicht möglich, mit Hilfe der tabellarischen Wertepaare eine numerische Abstandsminimierung von $\kappa(T_L)$ durchzuführen. Der Vorteil dieser Methode liegt in der genauen Abhängigkeit von $\kappa(T_L)$ im Entwicklungspunkt T₀ und in der anschließenden einfachen Integration.

Gleichung (3.1) ist die Differenzialgleichung für den Wärmestrom des Kristallgitters. Um eine für die Bauteilsimulation entsprechende Differenzengleichung angeben zu können, wird der Wärmestrom zwischen zwei Gitterpunkten x_i und x_j betrachtet, die die Temperaturen T_i bzw. T_j annehmen. Eine Näherung der Diskretisierung besteht darin, den Wärmestrom $\vec{S}_L(T_L)$ zwischen den Gitterpunkten als konstant anzunehmen. Dadurch ist es möglich, den Wärmestrom aus dem Integral herauszuziehen

$$\begin{eqnarray} \int_{x_i}^{x_j}\vec{S}_L(T_L)\cdot \mathrm{d} x=\vec{S}_L(T_i,... ...t \Delta x=-\int_{T_i}^{T_j}\kappa_L(T_L) \cdot \mathrm{d}T_L\; . \end{eqnarray}$$ (3.3)

Den Wert des Wärmestromes $\vec{S}_L$ zwischen zwei Diskretisierungspunkten erhält man mit (3.2) zu

$$\begin{eqnarray} \vec{S}_L(T_i,T_j)=-\frac{1}{\Delta x}\int_{T_i}^{T_j} C_1+C_2\cdot (T_L-T_0)+C_3 \cdot(T_L-T_0)^2\cdot \mathrm{d}T_L\; . \end{eqnarray}$$ (3.4)

Auswerten des Integrals und Zusammenfassen der Vorfaktoren zu den Potenzen liefert

$$\begin{eqnarray} \vec{S}_L(T_i,T_j)=-\frac{1}{\Delta x}\cdot (D_1\cdot(T_j-T_i) + D_2\cdot(T_j-T_i)^2 + D_3\cdot(T_j-T_i)^3) \end{eqnarray}$$ (3.5)

mit

$$\begin{eqnarray} D_1\!=\! C_1-C_2\!\cdot \! T_0+C_3\!\cdot\! T_0^2; \qquad D_2\!=\! C_2-2\!\cdot\! C_3\!\cdot \!T_0; \qquad D_3\!=\!C_3\; . \end{eqnarray}$$ (3.6)

Trägt man den berechneten Fluß in das Gleichungssytem ein, so ist er mit der entsprechenden Boxfläche A_i,j zu multiplizieren [18]. In MINIMOS-NT liegt die geometrische Information der Boxfläche bezogen auf den Punkteabstand in Form eines Flußfeldes F_i,j vor. Dabei beziehen sich die Indizes i,j auf die diskretisierten Gitterpunkte des Simulationsgitters. Sie besitzen einen definierten Abstand sowie eine zugeordnete Boxfläche. Bei einer zweidimensionalen Simulation entspricht die Fläche A_i,j einer Box-Seitenkante multipliziert mit der Einheitstiefe. Der Wärmefluß $\vec{S}_L$ vom Gitterpunkt i zum Gitterpunkt j berechnet sich daher zu

$$\begin{eqnarray} S_L^{i,j}=F_{i,j}\cdot \mathrm{M}_{\vec{S}_L}(T_i,T_j)\; . \end{eqnarray}$$ (3.7)

Die Größe $\mathrm{M}_{\vec{S}_L}$ ist die implementierte Modellfunktion des Wärmeflußes. Sie entspricht dem Term in der Klammer (rechte Seite von (3.5)) und berechnet neben der Information des Wärmeflußes auch die entsprechenden Ableitungen nach den Gittertemperaturen T_i,T_j.

Bei der Implementierung der Wärmeflußgleichung in MINIMOS-NT muß das verwendete Simulationsgitter berücksichtigt werden. Daher bezieht man alle Beiträge einer Gleichung auf das Boxvolumen, um entsprechend gewichtete Einträge zu erhalten. Die auftretende Divergenz einer Größe wird dabei mit Hilfe des GAUSS'schen Satzes (2.66) auf ein Oberflächenintegral zurückgeführt. In diskretisierter Form lautet die implementierte Energieflußgleichung im Drift-Diffusionsmodell für den Gitterpunkt i

$$\begin{eqnarray} \sum_j S_L^{i,j}=\sum_b \sum_j\frac{1}{2}\cdot\left[\frac{E_{b,... ...i-\psi_j)\right]\cdot I_b^{i,j}+V_i\cdot R_{net,i}\cdot E_{g}\; . \end{eqnarray}$$ (3.8)

Die Summation über j bedeutet eine Berechnung über alle angrenzenden Nachbarpunkte. Die Summation über b berechnet die Energieterme aller vorhandenen Ladungsträger. Die Größe R_i entspricht einer effektiven Rekombinationsrate. In MINIMOS-NT ist (3.8) durch eine sogenannte Kontrollfunktion f implementiert [18]

$$\begin{eqnarray} f_{\vec{S}_{L,DD,i}}=-\sum_j S_L^{i,j}+\sum_b \sum_j\frac{1}{2}... ...(\psi_i-\psi_j)\right]\cdot I_b^{i,j}+V_i\cdot R_i\cdot E_{g}\; . \end{eqnarray}$$ (3.9)

Der Faktor 1/2 ergibt sich durch die Zuordnung der Hälfte des Wärmeenergieeintrages zur Box i bzw. zur Box j. Der Fehler, der bei dieser Zuordnung gemacht wird, ist umso kleiner, je feiner das Simulationsgitter ist.

Die Kontrollfunktion für die hydrodynamische Gittererwärmung bestimmt sich laut (2.70) für den Gitterpunkt i zu

$$\begin{eqnarray} f_{\vec{S}_{L,HD,i}}=-\sum_j S_L^{i,j}+V_i\cdot\sum_b \frac{3\c... ...{T_{b,i} -T_{L,i}}{\tau_{w}}\right) -V_i\cdot\sum_b H_{eff,b}\; . \end{eqnarray}$$ (3.10)

Bei einer transienten Berechnung der Kristallgittererwärmung muß zusätzlich der zweite Term auf der linken Seite von (2.68) ausgewertet werden. Dabei wird eine Rückwärts-EULER-Zeitdiskretisierung, analog den elektrischen Gleichungen, angewendet [18]. Dies bedeutet, daß für den aktuellen Zeitschritt das Gleichungssytem mit den aktuellen Variablen berechnet wird. Nur für den transienten Term verwendet man die Variablen des letzten gerechneten Arbeitspunktes. Dadurch kann man den transienten Term von (2.68) folgendermaßen darstellen

$$\begin{eqnarray} \rho_L \cdot c_p \cdot \frac{\partial T_L}{\partial t}=\rho_L \cdot c_p(T_L,T_{L,old}) \cdot \frac{T_L - T_{L,old}}{\Delta t}\; . \end{eqnarray}$$ (3.11)

Im Gegensatz zur Massendichte $\rho_L$ ist die spezifische Wärmekapazität als c_p temperaturabhängig angenommen. Sie ist für die meisten Stoffe in den verlangten Temperaturbereichen relativ genau bekannt. Der Wert von c_p wird im Simulator durch eine Modellfunktion ausgewertet. Ähnlich wie bei der Wärmeleitung ist dafür im Simulator eine numerische Taylorreihe um $T_0\!=\!300\mathrm{K}$ implementiert

$$\begin{eqnarray} c_p(T_L,T_{L,old})=C_1+C_2\cdot T_{av}+C_3\cdot T_{av}^2 \quad \mathrm{mit}\quad T_{av} =\frac{1}{2}\cdot(T_L+T_{L,old})-T_0\; . \end{eqnarray}$$ (3.12)

Die transienten Kontrollfunktionen ergeben sich aus den stationären Kontrollfunktionen (3.9,3.10), jeweils erweitert um den transienten Term der Box i

$$\begin{eqnarray} V_i \cdot \rho_L\cdot c_p(T_L,T_{L,old}) \cdot \frac{T_L-T_{L,old}}{\Delta t}\; . \end{eqnarray}$$ (3.13)

Das Lösen des so aufgestellten Gleichungssystems berechnet Werte des neuen Zeitschrittes für die Gittertemperatur an den Gitterpunkten. Um eine Aussage über die Qualität der neuen Lösung treffen zu können, ist eine Abschätzung über den Fehler der Zeitdiskretisierung notwendig.

Bei rein thermischen Problemen, d.h. beim Lösen der Gleichung der Kristallgittererwärmung, kann man dafür die Änderung der Temperaturen der beiden Zeitschritte heranziehen. Man verlangt, daß der Abstand der Kristallgittertemperaturen zweier aufeinanderfolgender Zeitpunkte innerhalb einer festgelegten Schranke liegt. Dieser Abstand wird durch Normen der Lösungsvariablen (Kristallgittertertemperatur) bestimmt (siehe Kapitel 7.2.2). In welchem Bereich diese Abweichung liegen darf, hängt von der Genauigkeit der geforderten Lösung ab.

Löst man gleichzeitig zur Wärmeflußgleichung die Halbleitergleichungen, so hängt die Lösung der Variablen in den Halbleitergleichungen stark von der berechneten Gittertemperatur des entsprechenden Zeitschrittes ab. Beim Abschätzen der Qualität der transienten Simulation der Gittererwärmung kann man daher auch die Änderung des Potentials oder der Ladungsträgerkonzentrationen zweier aufeinanderfolgender Zeitschritte heranziehen. Dies setzt jedoch voraus, daß die angelegten Potentiale der betrachteten Zeitschritte konstant sind. Im Bauteilsimulator MINIMOS-NT kann man die thermische Gleichung mit den elektrischen Gleichungen gekoppelt oder ungekoppelt lösen. Bei der ungekoppelten Lösung löst man zuerst die Halbleitergleichungen bei festgehaltener Gittertemperatur. Anschließend hält man die elektrischen Größen fest und löst nach der Gittertemperatur. Dieser Vorgang wird so lange wiederholt, bis beim Wechseln der Gleichungssysteme die Abbruchbedingung erfüllt ist. Diese Lösungsmethode zeichnet sich durch ein stabiles Konvergenzverhalten aus, jedoch ist die Konvergenzgeschwindigkeit wesentlich langsamer als beim vollgekoppelten Lösen. Die entkoppelte Lösung stellt jedoch ein brauchbares Mittel dar, um das Konvergenzverhalten eines Gleichungssystems zu testen (siehe Kapitel 7). In MINIMOS-NT wird der vollgekoppelten Lösung der Vorzug gegeben, um das Gesamtsystem von thermischen und elektrischen Größen zu berechnen.