Stabilność algorytmu OFET Solver: błędy siatki, konwergencji i symulacji.

W tej sekcji wyjaśniono, jak diagnozować problemy ze stabilnością i zbieżnością rozwiązań w dwuwymiarowych symulacjach OFET. Symulacje jednowymiarowe są zazwyczaj stabilne i szybkie, ale wraz ze wzrostem wymiarowości z 1D do 2D i 3D liczba sprzężonych równań i powiązanych elementów gwałtownie wzrasta. Złe wybory siatki, nierealistyczne wartości ruchliwości, słabo zdefiniowane kontakty lub zbyt duże kroki polaryzacji mogą zatem powodować, że solver oparty na równaniach transportu ładunku przestanie generować zbieżne wyniki.

W tym miejscu celowo zastosujemy kilka prostych konfiguracji, aby wywołać sytuację, w której proces nie zbiega się do rozwiązania, a następnie zdiagnozujemy przyczyny tego zjawiska i sposoby jego naprawy. Odtwarzając samodzielnie te typowe scenariusze niepowodzeń, nauczysz się rozpoznawać objawy i stosować odpowiednie metody naprawcze (np. dostosowywanie siatki, zmniejszanie kroku korekcji lub modyfikowanie ekstremalnych parametrów materiału), aby uzyskać stabilne rozwiązanie.

1. Problemy ze stabilnością algorytmu spowodowane zbyt niską wartością współczynnika ruchliwości.

W wielu przypadkach problemy z konwergencją wynikają z użycia nierealnych wartości wejściowych. Aby to zilustrować, celowo ustawiamy ruchliwość nośników na wartość 0 i analizujemy wynikające z tego błędy. Oczywiście, w rzeczywistych materiałach ruchliwość nigdy nie jest dokładnie równa zeru – zawsze ma skończoną wartość – ale ten przykład pokazuje, co się dzieje, gdy algorytm otrzymuje nierealistyczny parametr. Ruchliwość jest konfigurowana w edytorze parametrów elektrycznych OFET. Jeśli wartość zostanie ustawiona poza fizycznie sensownym zakresem, równania dryfu i dyfuzji mogą stać się źle uwarunkowane lub algorytm może nie osiągnąć konwergencji.

Otwórz edytor parametrów elektrycznych (główne okno → zakładka Urządzenie) i ustaw wartość mobilności na 0, zgodnie z ilustracją ??. Następnie uruchom symulację. Powinien pojawić się ekran z komunikatem o błędzie, podobny do ??, z dwoma czerwonymi prostokątami: pionowym prostokątem otaczającym podsumowanie resztkowe oraz poziomym prostokątem w pobliżu komunikatu o błędzie.

Jeśli przyjrzysz się pionowemu, czerwonemu polu, zauważysz f()=2.54e2, a następnie inne linie, na których f() ma bardzo dużą wartość (np. 1e15, 1e5, 1e8, 1e7). W tym przypadku f() oznacza błąd resztkowy solvera: jest to miara tego, jak dobrze aktualnie spełnione są sprzężone równania. Idealnie, błąd resztkowy solvera powinien wynosić dokładnie zero, choć w praktyce nigdy nie udaje się tego osiągnąć. Wartość resztkowa w zakresie od 1e−10 do 1e−8 wskazuje na bardzo dobrą konwergencję, natomiast wartości w okolicach 1e−1 są nadal akceptowalne. Z kolei wartości resztkowe rzędu setek (np. 2.54e2) – zwłaszcza gdy towarzyszą im wartości komponentów tak duże jak 1e15 – są wyraźnym sygnałem, że solver ma problemy, a system stał się niestabilny.

Czerwony, poziomy pasek zazwyczaj wyświetla przyczynę, dla której symulacja zostaje przerwana: Holes asking for 3e5 but only defined in range [min, max]. Oznacza to, że obliczona wartość quasi-poziomu Fermiego wykracza poza zakres określony w tabeli. Przed rozpoczęciem symulacji model tworzy tabelę zależności, takich jak quasi-poziom Fermiego w funkcji gęstości nośników, w szerokim zakresie energii – zakres ten jest znacznie większy niż można by się spodziewać w przypadku typowego urządzenia. Gdy solver staje się niestabilny, model może zażądać nierealnych wartości, znacznie wykraczających poza ten zakres; w tym momencie solver nie może obliczyć wymaganych wartości i symulacja zostaje zatrzymana.

Krótko mówiąc, ustawiając wartość mobilności na nierealną wartość (μ = 0), zmusiliśmy równania do działania w reżimie, który nie ma sensu fizycznego, co powoduje, że metoda numeryczna zawodzi: pojawiają się duże błędy, a następnie program żąda wartości spoza wcześniej obliczonych tabel, co prowadzi do awarii. Ważne jest, aby zauważyć, że taki wynik nie jest wynikiem słabości modelu ani błędu w oprogramowaniu. Wprowadzając nierealne parametry, zmuszamy solver do działania w reżimie, w którym nie istnieje żadne matematycznie lub fizycznie poprawne rozwiązanie, dlatego nie może on dostarczyć sensownego wyniku.

2. Problemy ze stabilnością algorytmu spowodowane zbyt wysoką, nierealną wartością mobilności.

W tym przykładzie ustawiliśmy ruchliwość nośników na wartość μ = 1×106 m²·V⁻¹·s⁻¹. Jest to celowo nierealna wartość dla modelu dryfu i dyfuzji półprzewodnika: tak duża wartość może być bardziej związana z przewodnictwem metalicznym niż z transportem skokowym/przez pasy w organicznych lub typowych półprzewodnikach. W praktyce nie używalibyśmy tak dużej wartości w symulacjach OFET.

Jeśli uruchomisz symulację z tymi ustawieniami, możesz zaobserwować, że algorytm ma trudności z osiągnięciem zbieżności. W oknie wyświetlającym wartości reszt (okienko, które pokazuje f()), zobaczysz bardzo duże błędy – np. 1e24, 1e19, 1e21 – które nie maleją. To normalne, że wartości reszt są początkowo wysokie podczas pierwszych kilku iteracji, gdy algorytm szuka stabilnego stanu, ale powinny one szybko zmaleć do małych wartości. Uporczywie duże wartości reszt są wyraźnym sygnałem, że system stał się numerycznie niestabilny lub nie odpowiada rzeczywistości ze względu na wybrane parametry.

Wnioski: należy unikać nierealistycznie wysokich wartości mobilności. Zbyt duże współczynniki transportu powodują, że układ równań różniczkowych staje się bardzo wrażliwy, co prowadzi do słabego warunkowania i braku zbieżności. Należy stosować fizycznie uzasadnione wartości μ dla danego materiału i ponownie uruchomić symulację, stosując umiarkowane kroki zmiany parametrów, aby przywrócić stabilne zachowanie.

3. Nieprawidłowa konstrukcja elementów i problemy z połączeniami.

W tym przykładzie otwieramy edytor kontaktów i ustawiamy bardzo małą wartość szerokości kontaktu źródła (patrz ??), konkretnie jeden mikrometr. Na pierwszy rzut oka może się to wydawać nieistotne, ale powoduje to subtelny problem: symulacja opiera się na siatce różnic skończonych, w której wielkości takie jak pole elektryczne i gęstość nośników są obliczane tylko w zdefiniowanych punktach siatki.

Jeśli sprawdzisz edytor siatki elektrycznej (dostępny w zakładce „Elektryka”), zauważysz, że ta niezwykle cienka, jednomikronowa warstwa kontaktowa znajduje się pomiędzy punktami siatki (patrz ??). W rezultacie warstwa kontaktowa nie jest właściwie uwzględniana przez siatkę opartą na metodzie różnic skończonych i nigdy nie jest poprawnie aplikowana do struktury urządzenia. Powoduje to rozbieżność między zdefiniowaną geometrią a siatką numeryczną, co uniemożliwia osiągnięcie zbieżności symulacji.

Aby rozwiązać ten problem, można albo zwiększyć grubość obszaru styku, tak aby nakładał się on przynajmniej na jeden punkt siatki, albo zwiększyć gęstość siatki (czyli dodać więcej punktów) w siatce różnic skończonych, aby zapewnić prawidłowe odwzorowanie obszaru styku.

4. Gęstość siatki, zbieżność i nadmierne zagęszczenie siatki.