Προσομοίωση πολλαπλών χωρικών κλιμάκων και συστημική ανάλυση διεργασιών χημικής απόθεσης από ατμό

Abstract

Η διεργασία χημικής απόθεσης από ατμό (ΧΑΑ) αποτελεί μία από τις βασικές διεργασίες παραγωγής λεπτών στερεών υμενίων. Πραγματοποιείται σε αντιδραστήρες εξοπλισμένους με ειδικές επιφάνειες, τα δισκία (wafers), πάνω στις οποίες πραγματοποιούνται επιφανειακές χημικές αντιδράσεις και η απόθεση του υμενίου. Τα υμένια που παράγονται χρησιμοποιούνται από διατάξεις ημιαγωγών μέχρι μικρο- ή νανο-ηλεκτρομηχανικά συστήματα, νανο-ηλεκτρονικές διατάξεις και ολοκληρωμένα κυκλώματα. Η σύγχρονη μηχανική χαρακτηρίζεται από συνεχόμενη μείωση των διαστάσεων των κατασκευαζόμενων δομών και πλέον οι προδιαγραφές των παραγόμενων υμενίων αναφέρονται σε ιδιότητες μικρο- ή νανο-κλίμακας. Αν και η εξέλιξη των φυσικών/χημικών φαινομένων στις μικρο- και νανο-κλίμακες - που καθορίζουν τις ιδιότητες των υμενίων - εξαρτάται από τις (μακροσκοπικές) συνθήκες λειτουργίας του αντιδραστήρα, οι μίας-κλίμακας συμβατικές υπολογιστικές τεχνικές που έχουν αναπτυχθεί για την προσομοίωση διεργασιών ΧΑΑ, δεν μπορούν να “δουν” την εξέλιξη αυτή. Έτσι, καθίσταται απαραίτητη η προσομοίωση σε πολλαπλές χωρικές κλίμακες. Αναπτύσσονται δύο υπολογιστικά πλαίσια για την προσομοίωση πολλαπλών χωρικών κλιμάκων σε διεργασίες ΧΑΑ. Το πρώτο αφορά στη σύζευξη της μακρο-κλίμακας (τάξης cm) της αέριας φάσης του αντιδραστήρα ΧΑΑ με προσχηματισμένη μικρο-τοπογραφία (τάξης μm) στην επιφάνεια του δισκίου και το δεύτερο στη σύνδεση του κύριου όγκου του αντιδραστήρα ΧΑΑ με τη νανο-μορφολογία (τάξης nm) του αναπτυσσόμενου υμενίου. Η σύζευξη αφορά στην αμφίδρομη ενώ η σύνδεση στη μονόδρομη επικοινωνία μεταξύ των προτύπων στις διαφορετικές κλίμακες. Το πρότυπο (μοντέλο) της μακρο-κλίμακας είναι κοινό για τα δυο υπολογιστικά πλαίσια και περιλαμβάνει τις διαφορικές εξισώσεις με μερικές παραγώγους της διατήρησης μάζας, ορμής, ενέργειας και χημικών συστατικών εφόσον ισχύει η υπόθεση του συνεχούς μέσου (αριθμός Knudsen << 1). Η διακριτοποίηση και επίλυση των εξισώσεων για την ταχύτητα, πίεση, θερμοκρασία και συγκεντρώσεις συστατικών πραγματοποιείται με τη μέθοδο των πεπερασμένων όγκων ελέγχου και τον κώδικα υπολογιστικής ρευστοδυναμικής Αnsys 12/Fluent. To πρότυπο για τη μελέτη των φαινομένων σε προσχηματισμένη μικρο-τοπογραφία είναι ντετερμινιστικό και προκύπτει από τη σύζευξη τριών υπο-προτύπων: α) Βαλλιστικό πρότυπο υπολογισμού των τοπικών ροών των συστατικών μέσα στις δομές. β) Πρότυπο απόθεσης σε επιφάνειες, το οποίο εκφράζει την κινητική απόθεσης. γ) Αλγόριθμος εξέλιξης μετώπου βασισμένος στη μέθοδο των ισοϋψών. Ισχύει σε περιπτώσεις που η υπόθεση του συνεχούς μέσου καταρρέει (αριθμός Knudsen >> 1), όπως συμβαίνει μέσα στις δομές σε συνθήκες χαμηλής πίεσης. Το πρότυπο για την πρόβλεψη της νανο-μορφολογίας του αναπτυσσόμενου υμενίου είναι στοχαστικό και βασίζεται σε πρότυπο τύπου kinetic Monte Carlo (KMC). Η σύζευξη της μακρο-κλίμακας με τη μικρο-τοπογραφία στο δισκίο επιτυγχάνεται μέσω του συντελεστή ενεργού κατανάλωσης, ε, ο οποίος εισάγεται στην συνοριακή συνθήκη για την κατανάλωση των χημικών συστατικών στη αέρια φάση του αντιδραστήρα. Η ύπαρξη της μικρο-τοπογραφίας αυξάνει την ενεργό επιφάνεια στην οποία μπορεί να αποτεθεί υλικό με αποτέλεσμα την αύξηση της κατανάλωσης των χημικών συστατικών. Ο ε εισάγεται στη συνοριακή συνθήκη για τη διατήρηση των χημικών συστατικών προκειμένου η μακρο-κλίμακα να “αντιληθφεί” την αυξημένη αυτή κατανάλωση, και κατ’ επέκταση την ύπαρξη της μικρο-τοπογραφίας στο δισκίο, χωρίς η τελευταία να υπεισέρχεται στους υπολογισμούς της μακρο-κλίμακας. Κατά τη διάρκεια των υπολογισμών ο ε υπολογίζεται μέσω επαναληπτικής διαδικασίας. Η επικοινωνία μεταξύ των κωδίκων είναι αμφίδρομη και συνεχής καθόλη τη διάρκεια των υπολογισμών. Το υπολογιστικό πλαίσιο σύζευξης εφαρμόζεται σε διεργασία ΧΑΑ βολφραμίου (W) και πυριτίου (Si) με προσχηματισμένη μικρο-τοπογραφία αυλακιών στην επιφάνεια του δισκίου. Υπολογισμοί πραγματοποιούνται υπό διαφορετικές συνθήκες λειτουργίας του αντιδραστήρα προκειμένου να διερευνηθεί ο τρόπος πλήρωσης των αυλακιών. Η ανάλυση δείχνει ότι για υψηλές τιμές του συντελεστή προσκόλλησης των χημικών συστατικών που μετέχουν στις επιφανειακές χημικές αντιδράσεις, η απόθεση μέσα στα αυλάκια είναι ανισότροπη και προκαλεί τη δημιουργία κενού κατά την πλήρωση τους. Η ύπαρξη της μικρο-τοπογραφίας επηρεάζει το διάγραμμα Arrhenius της διεργασίας – δηλαδή το διάγραμμα που δείχνει την εξάρτηση του ρυθμού απόθεσης του υμενίου από τη θερμοκρασία και οριοθετεί τις περιοχές: Περιοχή αντίδρασης, όπου ελέγχων μηχανισμός του ρυθμού απόθεσης είναι η χημική αντίδραση, περιοχή φαινομένων μεταφοράς ή διάχυσης, όπου ελέγχων μηχανισμός του ρυθμού απόθεσης είναι η διάχυση των χημικών συστατικών προς το δισκίο και μεταβατική περιοχή στην οποία και οι δύο παραπάνω μηχανισμοί είναι σημαντικοί. Συγκεκριμένα, η ύπαρξη της μικρο-τοπογραφίας αυξάνει την ενεργό επιφάνεια για απόθεση αυξάνοντας την κατανάλωση των αντιδρώντων αλλά προκαλώντας τη μείωση του ρυθμού απόθεσης. Η μείωση οφείλεται στο φαινόμενο εξάντλησης (loading phenomenon), το οποίο είναι ιδιαίτερα σημαντικό στη μεταβατική περιοχή και στην περιοχή διάχυσης. Έτσι, εξαιτίας της ύπαρξης της μικρο-τοπογραφίας και μέσω του φαινομένου εξάντλησης, τα όρια της μεταβατικής περιοχής και της περιοχής διάχυσης στο διάγραμμα Arrhenius μετακινούνται σε χαμηλότερες θερμοκρασίες. Επιπλέον, η ύπαρξη της μικρο-τοπογραφίας στο δισκίο επηρεάζει την αέρια φάση του αντιδραστήρα, λόγω της αυξημένης κατανάλωσης των αντιδρώντων, μεταβάλλοντας την κατανομή των κλασμάτων μάζας μέχρι και την είσοδο του αντιδραστήρα. Η σύνδεση της μακρο-κλίμακας με τη νανο-μορφολογία του αναπτυσσόμενου υμενίου επιτυγχάνεται θεωρώντας ότι ο ρυθμός απόθεσης είναι ανεξάρτητος από την κλίμακα στην οποία υπολογίζεται. Έτσι, δεν είναι αναγκαία η συνεχής ανταλλαγή πληροφορίας μεταξύ των κωδίκων: Μόνο ο κώδικας της μακρο-κλίμακας τροφοδοτεί με πληροφορία - τον ρυθμό απόθεσης - τον κώδικα KΜC για την πρόβλεψη της νανο-μορφολογίας του αναπτυσσόμενου υμενίου. Η μεθοδολογία σύνδεσης εφαρμόζεται σε διεργασία ΧΑΑ Si λαμβάνοντας υπόψη την ύπαρξη των διμερών, χαρακτηριστικών δομών στην επιφάνεια Si. Οι υπολογισμοί δείχνουν ότι κατά μήκος της ακτίνας του δισκίου μπορεί να παρατηρηθεί διαφορετικός προσανατολισμός των διμερών εξαιτίας του διαφορετικού ρυθμού απόθεσης. Η περιγραφή και ανάλυση φυσικών/χημικών φαινομένων στις πολλαπλές κλίμακες δεν απαιτεί μόνο σύνθετες “πολυδιάστατες” προσεγγίσεις αλλά και αποτελεσματικές παράλληλες μεθοδολογίες προκειμένου να καλυφθούν οι υπολογιστικές απαιτήσεις των προσεγγίσεων αυτών τόσο σε επίπεδο υπολογιστικού χρόνου όσο και σε μνήμη. Υλοποιείται μία υβριδική πολυ-παράλληλη μέθοδος, η οποία συνδυάζει αφενός μεθόδους διαμοιρασμού χωρίου για τους υπολογισμούς στη μακρο-κλίμακα του αντιδραστήρα και αφετέρου την παράλληλη μέθοδο “αφέντη-εργάτη” για τους υπολογισμούς με το πρότυπο της μικρο-κλίμακας. Ο όρος “πολυ-παράλληλη” προκύπτει από το ότι διαφορετικός αριθμός επεξεργαστών μπορεί να χρησιμοποιηθεί για την επίλυση των προβλημάτων στις διαφορετικές κλίμακες. Το τελευταίο είναι ιδιαίτερα σημαντικό καθώς οι υπολογιστικές απαιτήσεις στις διαφορετικές κλίμακες είναι διαφορετικές. Χρησιμοποιώντας μόνο τη μέθοδο αφέντη-εργάτη μπορούν να επιτευχθούν επιταχύνσεις κοντά στην ιδανική μειώνοντας αποτελεσματικά τον υπολογιστικό χρόνο, που σε ορισμένες περιπτώσεις από τις 3 μέρες μειώθηκε στις 3 ½ ώρες. Επιπλέον, υπολογισμοί με την υβριδική πολυ-παράλληλη μέθοδο δείχνουν ότι οι διαφορετικές παράλληλοι μέθοδοι “συνεργάζονται” αποτελεσματικά μεταξύ τους προκειμένου να επιταχυνθούν οι υπολογισμοί σε όλο το εύρος των κλιμάκων. Για την αριθμητική επίλυση των εξισώσεων που διέπουν προβλήματα φαινομένων μεταφοράς, η χρησιμοποίηση ώριμων εμπορικών κωδίκων - όπως ο Fluent - έχει γίνει κοινή πρακτική. Παρόλη την ανάπτυξή τους, οι εμπορικοί αυτοί κώδικες, δεν είναι ικανοί να πραγματοποιήσουν επαρκή συστημική ανάλυση, δηλαδή συστηματική αναζήτηση και εντοπισμό πολλαπλών λύσεων και σημείων στροφής πάνω σε κλάδους λύσεων και να υπολογίσουν ασταθείς λύσεις μόνιμης κατάστασης συναρτήσει των μεταβαλλόμενων παραμέτρων. Έτσι, αδυνατούν να προσφέρουν “όλα τα κομμάτια από το πάζλ”, δηλαδή την εξάρτηση των λύσεων μη γραμμικών προβλημάτων από κρίσιμες παραμέτρους. Χαμένα κομμάτια μπορεί να κρύβουν σημαντική πληροφορία για τα όρια της ευστάθειας της λύσης όπως επίσης και ολόκληρους κλάδους οι οποίοι μπορεί να είναι επωφελείς για μια διεργασία. Αναπτύσσεται υπολογιστικό πλαίσιο που βασίζεται στη μέθοδο αναδρομικής προβολής (recursive projection method ή RPM) και εφαρμόζεται σαν εξωτερικό υπολογιστικό κέλυφος γύρω από τον Fluent. Ο κύριος σκοπός είναι να εξαναγκάσει τον Fluent να συγκλίνει σε κλάδους λύσεων μη γραμμικών προβλημάτων με συστηματικό και αποτελεσματικό τρόπο ακόμη και σε ασταθείς μόνιμες καταστάσεις. H μέθοδος RPM/Fluent χρησιμοποιείται με επιτυχία για να βρει έναν ολόκληρο χώρο λύσεων σε διεργασία ΧΑΑ Si. Ο χώρος λύσεων αποτελείται από έναν ευσταθή κλάδο, στον οποίο κυριαρχεί η φυσική συναγωγή, ακολουθούμενος, μετά από ένα πρώτο σημείο στροφής, από έναν ασταθή κλάδο ο οποίος οδηγεί, έπειτα από ένα δεύτερο σημείο στροφής, σε ένα ευσταθή κλάδο λύσεων στον οποίο ο κυρίαρχος μηχανισμός συναγωγής είναι η εξαναγκασμένη συναγωγή. Χαρακτηριστικοί ρυθμοί απόθεσης παρουσιάζονται σε κάθε κλάδο λύσεων κατά μήκος το δισκίου και υπολογίζεται η ανομοιομορφία τους. Τα αποτελέσματα δείχνουν ότι είναι σημαντική η γνώση της πολλαπλότητας λύσεων για τις ίδιες τιμές παραμέτρων λειτουργίας καθώς το τελικό προϊόν, το υμένιο, διαφέρει ανάλογα με τον κυρίαρχο μηχανισμό συναγωγής. Για τιμές παραμέτρων στις οποίες συνυπάρχουν πολλαπλά ροϊκά πεδία, υπάρχει σημαντική διαφορά στους ρυθμούς απόθεσης. Δύο διαφορετικά διαγράμματα Arrhenius υπολογίζονται για τις ίδιες συνθήκες λειτουργίας του αντιδραστήρα που αντιστοιχούν στη φυσική και στην εξαναγκασμένη συναγωγή.

The produced films via chemical vapor deposition (CVD) are utilized to a wide range of applications; from semiconductor devices to micro- and nano-electromechanical systems (MEMS and NEMS), and protein microarrays and chips. Nowadays, the size of these devices or systems shrinks to lower scales and the specifications of the films, e.g. thickness, conformality (thickness uniformity on a patterned surface), surface morphology, refer to properties in micro- or nano-scale. Thus, the single scale conventional CVD modeling methods are not adequate and more advanced, multiscale modeling, methods are needed for studying the phenomena in the co-existing (multiple length) scales. For example, the filling of a micro-trench on the wafer and the pertinent film conformality come from the “interaction” of the macro- or reactor scale with the micro- or feature (trench) scale. Similarly, the nano-roughness developing on a film’s surface during deposition comes from the interaction of the macro- or reactor scale with the submicro- or nano- or surface morphology scale. The description of each scale requires a model: The reactor scale model (RSM) is used for the description of the transport phenomena in the bulk phase of the CVD reactor, the feature scale model (FSM) is used to describe the film deposition inside the features (e.g. trenches) and the nano-morphology model (NMM) is used to trail the surface morphology of the deposited film. The effect of the operational parameters of a CVD reactor on the film profile evolution inside a trench on the wafer or on the film’s surface nano-morphology evolution can be predicted by linking or coupling of RSM with FSM or NMM. The linking of models refers to their sequential use, while the coupling incorporates a two-way interaction of the models. A multiscale computational framework is described and applied to CVD of a film on a predefined topography, i.e. an array of trenches on the wafer; the latter computations require the coupling of RSM with FSM. Additionally, the framework is applied to CVD of a film on an initially flat wafer; linking of the RSM with NMM is performed to calculate the surface morphology of the deposited film. The multiscale computational framework consists of a RSM, a FSM, a NMM and the algorithms for coupling RSM with FSM and linking RSM with NMM. The RSM describes the transport phenomena in the macro-scale of the CVD reactor. The governing equations are the continuity, the momentum, the energy and the species transport equations. They are solved numerically at steady state to predict the velocity, pressure, temperature and species concentrations inside the bulk phase of the CVD reactor. The CFD code Ansys 12/Fluent is used for the numerical solution of the aforementioned set of equations. The FSM results from coupling a local flux calculation model, a surface model and a profile evolution algorithm. The local flux calculation model is a ballistic model, which is formulated by a set of integral equations and it is valid at the high Knudsen number conditions occurring in the micro-trenches of the predefined topography. It links the fluxes of the species arriving just above the predefined topography on the wafer with the local fluxes inside the features (e.g. trenches) of the topography. The surface model describes the surface processes and quantifies the effect of species’ fluxes or concentrations on the local deposition rate. The profile evolution algorithm of the deposited film inside the features (moving boundary) uses the level set method. The NMM is based on the kinetic Monte Carlo method. All surface processes, i.e. adsorption, surface diffusion, surface reaction, desorption are modeled as a Markov process by transition probabilities per unit time. The transition probabilities per unit time are modeled by an Arrhenius expression and depend on local activation energies. The NMM is essentially a surface model as the surface model of the FSM; however, the surface model of the FSM is formulated by an analytical surface reaction rate expression (e.g. in terms of the reactive species concentration) and cannot predict the surface morphology. The NMM can also “operate” at sub-micro- or even at micro-scale by the implementation of specific coarse grained methods. The coupling of the RSM with FSM is based on the correction of consumption rates of each species on a predefined topography (e.g. an array of trenches) on the wafer. The aim of this correction is to take into account the increased consumption of species inside the topography, without the topography being included in the computational domain of macro-scale. A correction factor, εk, is applied to each surface reaction rate k, reflecting a change of the boundary condition of the species equation. The convergence of the iterative scheme essentially ensures a mass conservation for the species. The linking of the RSM with NMM starts with the numerical solution of the equations in the macro-scale at steady state. All species fluxes just above the wafer are fed to NMM. No bidirectional exchange of computational information between the scales is performed. In other words, no effect of the nano-scale on the macro-scale is considered. This is a reasonable assumption given that the change of the surface nano-morphology is not expected to alter the species consumption on the wafer surface. To describe the physical and chemical phenomena in the macro-scale of the CVD reactor, the RSM based on the partial differential equations (PDEs) of the conservation of mass, momentum, energy and chemical species, is used. For numerically solving the equations, Fluent solver is used along with the appropriate boundary conditions (BCs). The appropriate BCs for the aforementioned set of PDEs come from the solution of a problem in a different co-existing scale, i.e. the micro- or nano-scale on the wafer’s surface. The description of such complex phenomena demands not only sophisticated multidisciplinary approaches but also efficient parallel methodologies that can handle the associated high computational demands. We propose a crossbred multi-parallel method (CMPM) by combining (crossbreeding) different parallel techniques in the multiple scales of interest. CMPM is demonstrated on the coupling methodology, described in the previous section, for coupling a RSM with FSM. In particular, we combine domain decomposition techniques to accelerate the computations with RSM which is based on PDEs, along with a “master-worker” scheme to calculate the appropriate BCs for the aforementioned PDEs. These BCs come from computations with FSM and are changing during the simulation time. The use of master-worker parallel scheme is feasible since the computations with FSM can be performed independently for each boundary cell of the wafer where a cluster of trenches is adjacent [see Fig. 5(b)]. The multi-parallel term stems from the fact that different number of processors can be used to implement the domain decomposition and the master-worker schemes, an important factor since the computational demands in the different scales varies. The applicability of CMPM is pronounced in problems where the computations in the different, multiple scales are performed independently from each other. The potential of applying the CMPM in problems where special treatment of the BCs is needed, such as in fluid-structure interaction problems, is under investigation. Nowadays, computational fluid dynamics (CFD) codes are becoming standard in many fields of science and engineering involving flows of gases and liquids; numerical simulations are used from plasma processing and electromagnetics to catalysis and chemical vapor deposition processes. Despite their evolution, commercial CFD codes are practically unable to perform adequate systemic analysis, that is systematic detection and tracing of multiple solutions, to circumvent turning point singularities along solution branches and to compute unstable steady-state solutions as parameter vary. Therefore, they fail to provide all the “pieces of the puzzle”. i.e. of the dependence of solutions of nonlinear problems on key parameters. These shortcomings are serious limitations on the predictive capability of commercial CFD codes: missing pieces may hide crucial information about the limits of stability of solutions as well as entire branches of solutions which might be suggestive of advantageous operating “windows” of the process of concern. We propose a framework for efficiently enabling a CFD code, in this case Ansys 12/Fluent to converge past turning points on unstable steady states branches and therefore successfully trace an entire solution branch. The proposed framework is based on the so-called recursive projection method (RPM) which is implemented as a computational shell around the CFD commercial code Fluent. The main purpose is to enable Fluent to trace solution branches of nonlinear problems that have multiple solutions in a systematic and efficient way, and even induce convergence on unstable steady states. The RPM is combined with arc-length type continuation techniques in order to enable convergence on unstable solutions past turning points, offering therefore the possibility to trace entire solution branches. As a by-product it also delivers approximations of the critical eigenmodes of the discretized physical problem, through the solution of a low-dimensional eigenvalue problem. This proves useful in problems with bifurcation points, where we can then use the derived eigenvectors for solution branch switching. Finally, it helps reduce computational cost, when used in the context of parameter continuation.