Μαθηματική προσομοίωση και μελέτη μαγνητοϋδροδυναμικών ροών σε καμπυλόγραμμες γεωμετρίες

Abstract

The main object of the present thesis is the study of magnetohydrodynamic (MHD) flows with inductionless approaches, using numerical models of computational fluid dynamics. Magnetohydrodynamics deals with flows of electrically conducting fluids which are subject to a magnetic field and/or an electric current driven by an external voltage. Thus, in the fluid there is a dynamically and mutually interaction of inertial, viscous and electromagnetic forces under the presence of a magnetic field. The research is focused in regions where a strong coupling of fluid mechanics and electromagnetism is observed. These regions are characterized by a mighty non-linearity and are of great importance in the designing of equipment for fusion reactors concerning especially the cooling devices. Thus, one of the main purposes of this study was the development of capable numerical models in order to predict the MHD behaviour of the conducting fluids with the lowest possible computational cost. Another motive was the proper implementation of them in areas, where the asymptotic-analytic tools are not valid (i.e. low Ha numbers), in order to yield reliable results after having been tested and validated in ranges of high and moderate Ha numbers where it is possible. The particularities of computationally handling MHD flows are, the enormous source terms, that increase exponentially with the Ha parameter and their heterogeneity throughout the computational domain (high Lorentz forces appears in very thin Hartmann and side layers or in tangential layers, while simultaneously there are practically absent in the core flow). Concurrently, the small convection or/and diffusion terms, which appear for high Ha, imbalance the discretised equations having as a result the attenuation of the main diagonal in the solution matrix of linear system. In addition, the many and intensive discontinuities of the computational domains and the large non-uniform grids impose the use of low relaxation factors expanding the total execution time. Based on the preceding, four semi-implicit algorithms are evaluated for pressure-velocity coupling of the MHD transport equations. Form the comparative study follows that the scheme with the faster convergence response in such problems, where coexist large source terms and rigorous dissimilarities (like, natural convection and MHD natural convection flows), is the SIMPLEC algorithm, while the best robustness is presented by the SIMPLEX. Next important issue is the solution methods of linear systems deriving from elliptic PDEs, comprising the main body of numerical simulation software. Normally, the program spends most of its execution time in solving the linearized transport equations, and so the efficiency of the linear solvers underpins completely the efficiency of the solution method. Therefore, a crucial aspect of any efficient solution of such a problem is the speed that these linearized equations are solved taking account also the MHD singularities. The choice of a fast capable method assists to reduce the time execution, importantly. A number of different algorithms, originated by four algorithm families (direct & iterative methods, static, ADI, strongly implicit and Krylov space methods) are contrasted and their resource use is compared in terms of speed and secondly of memory use. Good and capable characteristics appear the family of strongly implicit solution methods (based on incomplete LU factorization), while the Conjugate Gradient methods present excellent efficiency. The convergence's rate of Krylov space methods can be improved by the use of a preconditioner, which can be probably a method of the strongly implicit family. Simultaneously, the discretisation schemes of the convection terms and temporal variations are investigated and evaluated by the use of kwon benchmark tests prejudging third order schemes, which belong to TVD family schemes. These schemes preserve the property to reduce their accuracy to lower order in areas with discontinuities in order to maintain their monotonicity. Oppositely, this cannot be done by the common schemes like QUICK, SOU and CD. The combination of space and temporal discretisation schemes is also experienced in order to be drawn conclusions for the selection of the proper combination limiting to the minimum the numerical errors. Additionally, some specific and appropriate modifications which might be considered in a CFD code for a MHD simulation are also discussed. The first of them is the discretisation of the Lorentz force term. From the numerical point of view, because the Lorentz force is greatly increased with increasing Hartman number, and since it can be larger than other terms in the governing equation, its treatment may control the accuracy and convergence of the MHD calculation. However, most calculation efforts have concentrated on discretization of the convection term and the pressure-velocity coupling algorithm, even in the high Hartman number MHD flow. However, the existence of an analytic solution of a simple MHD flow, which is independent of the imposed magnetic field, helps to be proven that the Lorentz cannot be handled as a simple source term because the common treatment leads to inaccurate results. The numerical error is analogous to the square of Ha parameter. Although, the method that produces the most accurate approach is computationally expensive and demanding. Another interesting point, which is related to the electric potential equation, is the calculation of it in 3D domain with insulating or finite conducting walls. The solution of the electric potential demands the additional resolving of a Poison equation on every boundary surface in order to be satisfied the current conservation in the computational domain. This requires an extra iterative procedure in order to equate the normal and the tangential currents on every boundary surface presenting the latter electric properties. Simulations of MHD natural convection compose a preferential area for the proper validation of the computational methods due to numerous studies and papers related to the above problem and techniques. Thus, the successful comparison of the upper tools with the published data creates a reliable base for further explorations to other areas of computational MHD. One of these is the study of a MHD rotating flow enclosed in cylindrical and spherical shells. This analysis considers two basic flow configurations covering a variety of investigations described in literature. In the first case, a fluid is enclosed in a cylinder with electrically perfectly conducting walls and the flow is driven by a small rotating, conducting disk. In the second case, a flow is considered in a spherical shell with an inner rotating sphere. The fluid in both cases is subjected to an external axial uniform magnetic field. The results show that these flows exhibit two different types of cores, separated from each other by a tangential layer parallel to the axis of rotation. The inner core follows a solid-body rotation while the outer is quasi stagnant. A counter-rotating jet is developed which is located in the tangential layer between the cores. The characteristics of the tangential layer are determined and the properties of the secondary motion are investigated for a wide range of Hartmann numbers. The main aim of present work is to enlarge the existing knowledge about such flows, and to support basic assumptions made in the past for treating this kind of flows by asymptotic techniques. The main features which had been discovered in the past for differentially rotating shells are recovered. In particular, the properties of the counter rotating jet that forms at the so-called tangent cylinder (parallel to magnetic field lines) at the edge of the rotating body are investigated in detail. The dependence on the strength of the magnetic field, i.e. the Hartmann number, of the jet velocity and thickness, and of the intensity of the secondary motion driven by centrifugal forces, is assessed. For high Hartmann numbers, the present results exhibit two different types of flow cores. One is located inside the tangent cylinder that is oriented parallel to the magnetic field lines between the rotating body and the outer wall, and the other outside it. Both cores are separated from each other by a thin, viscous tangential layer, in which the fluid reaches highest circumferential velocities. The inner core follows a solid body-type rotation with the fluid near the perfectly conducting body rotating at the same angular velocity. Its velocity decays linearly in the axial direction and vanishes at the outer non-moving wall. The outer core instead is practically stagnant. In the case of the rotating disk, the Hartmann layers are absent since the obtained core solution satisfies already the no-slip boundary conditions at the disk and the outer walls. In the case of the rotating sphere, instead, the inner core near the axis rotates faster than the sphere and slower closer to the tangent cylinder. This is caused by the curvature of the inner and outer boundaries and by the Hartmann length (height of the gap measured along magnetic field lines), which varies with radius r. As a result Hartmann layers appear at the rotating sphere, across which the no-slip condition is satisfied. For both cases (disk and sphere), the raising of Ha reduces the thickness of the tangential layer as δ~Ha⁻¹/² and increases its angular velocity of counter-rotation like ω~Ha¹/². The kinetic energy of the secondary motion driven by centrifugal forces is reduced proportional to Ha parameter. At the same time, the electric potential develops strong gradients in the transitional "column" that is formed at the tangent cylinder. The distribution of potential in the inner core follows a linear decay in the axial direction and vanishes at the external boundary. The driving mechanism of the rotating flow is the radial potential gradient that is observed in the inner core. However, the highest potential gradient (of opposite sign to that in the core) is observed in the symmetry plane at the outer periphery of the rotating body, there causing the highest counter-rotating velocities. For the high Hartmann numbers considered, all quantities calculated for the cores are in good agreement with results obtained by asymptotic theory. The results obtained in the tangential layer differ slightly from those derived by analytic methods for the case of the rotating disk but they confirm well the asymptotic behaviour. Very good agreement is achieved for the rotating sphere for which a proper numerical resolution of the Hartmann layer requiring 10-12 cells across the layer thickness could be achieved easier, since this layer is not as thin as in the disk, due to the curvature of the sphere. Although, in this kind of flows the Hartmann layers do not carry a significant amount of current, their proper numerical resolution turns out to be the key for a successful and accurate simulation. The present work complements known results for perfectly conducting walls, and to support the asymptotic theory, usually applied for high Hartmann numbers, where secondary flows are neglected a priori. For lower Hartmann numbers, however, where the asymptotic theory is not valid, the present numerical simulations yield reliable results. As second main research topic is considered the study of the walls conductivity effect on the flow of a conducting fluid sheltered in a cylindrical cavity. In this configuration the MHD flow is driven by a rotating disk at the top of a cylinder. The effect of an axial magnetic field on the flow is investigated for an aspect ratio H/R equal to 1. The magnetic Reynolds number is assumed to be small whereas the interaction parameter, N, is large compared to unity. This allows deriving asymptotic results for the flow solution. Various combinations of the top, bottom and vertical walls conductivities are simulated. The results obtained showed that one can control the primary flow through a good choice of the electrical conductivity of both the disk and cylinder walls. Also, in the presence of a strong vertical magnetic field, the computed flow solution, electric potential and current distribution have been compared to asymptotic solutions which can be obtained when inertial effects are neglected and the agreement is good. The influence of the electrical conductivity of the walls has been also considered. The results indicate that one can control the behaviour of the primary flow through a good choice of the electric conductivities of the disk and cylinder walls. In view of possible future laboratory experiments aiming at investigating the effect of a strong magnetic field on chemical phenomena occurring at the interface between a solid wall and a liquid metal flow the simple geometry studied here appears to be a convenient tool which offers the possibility to vary velocity gradients by modifying the wall conductivities. In another study, the effects of the magnetic field, the fluid and wall electrical conductivities, and the wall thickness are investigated. The study was performed for constant Re=100 and for Ha in the range 0≤Ha≤100. This corresponds to a range of interaction parameter N of 0≤N≤100. The work focuses on thin walls, which simplifies the boundary conditions. The thin wall boundary condition is used for the first time for a moving wall. It is shown that for fixed values of the Hartmann and Reynolds numbers, the velocity distribution depends strongly on the conductance ratio k, in spite of the fact that, the Hartmann layer thickness and side layer thickness do not vary with k. The study is destined to predict the influence of a magnetic field on the corrosion rate of a liquid metal on a metallic wall. It is assumed that this corrosion is controlled by the near-wall hydrodynamic which is then controlled by an external magnetic field. The concentration equation for the corrosion product is solved, and predicts the evolution of the mass transfer with Ha. At same magnitude of Ha the mass transfer is higher for conducting than insulating walls. The calculations for large values of the conductance ratio k indicate that the hydrodynamics is largely controlled by this factor. Even for very small values of k, the wall cannot be assimilated to be insulating case when the Hartmann number is high enough when, even for k < 1, the flow is very close to the perfectly conducting case. This can be easily explained by the sharing of the electric current between the walls, and the Hartmann and parallel layers. Depending on the thickness of the wall (taken into account by the conductivity ratio) and the thickness of the boundary layers, the electric current closes preferentially inside the wall or inside the layers. For large values of Ha, the two layers are very small and then the current closes preferentially in the walls even if they are very thin. The azimuthally velocity is classically organized into core, Hartmann layers, and a parallel layer. The core flow exhibits a linear velocity variation in the axial direction, with a slope that depends strongly on k. The boundary layers thickness is constant with k, but the azimuthally velocity gradient at the rotating disk (and thus the mass transfer process) varies considerably with this factor. The electric current lines are parallel to the imposed magnetic field in the core flow, and at the vicinity of both disks for the perfectly conducting case. Thus, the vanishing of the Lorentz force explains why the meridian flow rate is stronger for high k. Numerical simulations of mass transfer, for high Schmidt number, from the rotating disk in laminar flow show that magnetic field changes significantly the phenomena. The presence of the magnetic field in the present geometry damps the mass transfer. The electric properties of the walls, through the conductance ratio k, play a significant role in the mass transfer. So, if the corrosion processes is controlled by the near walls hydrodynamics, the conductance ratio of the walls is an important factor which cannot be neglected.

Αντικείμενο της παρούσας διατριβής είναι η μελέτη μαγνητοϋδροδυναμικών ροών με τη χρήση μαθηματικών μοντέλων υπολογιστικής ρευστοδυναμικής υπό το πρίσμα μη επαγωγικών μεθόδων. Ως μελέτη ΜΥΔ ροών ορίζεται η μελέτη των αμοιβαίων αλληλεπιδράσεων μεταξύ των δυνάμεων αδρανείας, ιξώδους και ηλεκτρομαγνητικών δυνάμεων σε στοιχεία αγώγιμου ρευστού εκτεθειμένα σε μαγνητικά πεδία. Σκοπός της ΜΥΔ ανάλυσης είναι να αναδείξει τον συσχετισμό δυνάμενων που συνυπάρχουν κατά τη ροή αγώγιμου ρευστού σε συγκεκριμένες δυναμικές περιοχές. Η μελέτη επικεντρώνεται στην περιοχή στην οποία εμφανίζεται ισχυρή σύζευξη ρευστομηχανικών και ηλεκτρομαγνητικών μεγεθών. Περιοχή η οποία, αφενός μεν χαρακτηρίζεται από ισχυρή μη-γραμμικότητα αφετέρου δε, είναι ιδιαίτερης σημασίας για το σχεδίασμά εξαρτημάτων των υποσυστημάτων ψύξης για αντιδραστήρες σύντηξης. Βασικός στόχος της παρούσας διατριβής είναι η δημιουργία μαθηματικών μοντέλων ικανών να προβλέπουν την ΜΥΔ συμπεριφορά αγώγιμων ρευστών με το μικρότερο δυνατό υπολογιστικό κόστος και η εφαρμογή τους σε περιοχές όπου οι ασυμπτωτικές-αναλυτικές μέθοδοι δεν ισχύουν. Οι ιδιαιτερότητες κατά την προσομοίωση ΜΥΔ ροών είναι, οι μεγάλοι όροι πηγής που αυξάνουν εκθετικά με την αύξηση του μαγνητικής επαγωγής και η ανομοιομορφία των συγκεκριμένων όρων στο εύρος τιμών τους (εμφάνιση υψηλών τιμών σε λεπτόπαχα Hartmann και πλάγια στρώματα ή σε επίσης λεπτές παράλληλες δέσμες, με ταυτόχρονα μικρές τιμές των ίδιων όρων στους πυρήνες της ροής). Ταυτόχρονα, οι αρκετά μικροί όροι συναγωγής ή και διάχυσης, που δημιουργούν ανισορροπία κατά την επίλυση των διακριτοποιημένων εξισώσεων, έχουν ως αποτέλεσμα την εξασθένιση της κυρίας διαγώνιου στο διαμορφούμενο αραιό μητρώο επίλυσης των γραμμικοποιημένων εξισώσεων. Παράλληλα, οι πολλές και έντονες ασυνέχειες στο πεδίου υπολογισμών και τα μεγάλα υπολογιστικά πλέγματα επιβάλλουν μικρούς συντελεστές υποχαλάρωσης που οδηγούν στην επιμήκυνση του χρόνου εκτέλεσης. Συνεπώς με βάση τα παραπάνω, αξιολογούνται αρχικά, τέσσερις βασικοί αλγόριθμοι ημι-πεπλεγμένης σύζευξης των βασικών εξισώσεων ΜΥΔ ροών. Από την συγκριτική δοκιμασία σύζευξης πίεσης-ορμής προκύπτει ότι το σχήμα με την ταχύτερη απόκριση για την επίτευξη σύγκλισης σε προβλήματα ΥΡΔ με μεγάλους όρους πηγής και σε πεδία που εμφανίζουν έντονη ανομοιογένεια (π.χ. προβλήματα φυσικής συναγωγής και ΜΥΔ) είναι ο αλγόριθμος SIMPLEC, ενώ την καλύτερη προσαρμοστικότητα εμφανίζει ο αλγόριθμος SIMPLEX. Στην συνέχεια αναλύονται και συγκρίνονται οι αλγόριθμοι επίλυσης γραμμικών συστημάτων που αποτελούν τον βασικό κορμό ενός λογισμικού αριθμητικής προσομοίωσης πεπλεγμένων εξισώσεων. Η ανάγκη για την ταχύτερη επίλυση μεγάλων συστημάτων γραμμικών εξισώσεων που προκύπτουν από την διακριτοποίηση των ΜΔΕ είναι επιτακτική και εάν σε αυτό προστεθούν, και οι ιδιαιτερότητες της ΜΥΔ, τότε η διερεύνηση αυτών των μεθόδων αποκτά ιδιαίτερη σημασία με αποτέλεσμα η επιλογή ενός αποδοτικού και γρήγορου αλγόριθμου βοηθά σημαντικά στην περιστολή του χρόνου εκτέλεσης. Για τους παραπάνω λόγους αντιπαραβάλλονται τέσσερις οικογένειες αλγορίθμων επίλυσης γραμμικών συστημάτων (στατικές, κατευθυντικής αποικοδόμησης, ισχυρά πεπλεγμένες μέθοδοι και μέθοδοι συζυγών κλίσεων). Βέλτιστα χαρακτηριστικά εμφανίζουν οι αλγόριθμοι της ομάδας των ισχυρά πεπλεγμένων μεθόδων και άριστη συμπεριφορά η ομάδα των συζυγών κλίσεων. Οι επιδόσεις της τελευταίας ομάδας μπορούν να βελτιωθούν με την χρήση προρυθμιστών και να αποτελέσουν αποδοτικό εργαλείο στην επίλυση μεγάλων γραμμικών συστημάτων. Συγχρόνως, διερευνώνται σχήματα διακριτοποίησης των όρων συναγωγής και χρονικής μεταβολής. Η αξιολόγηση των σχημάτων διακριτοποίησης γίνεται με γνωστές από την βιβλιογραφία δοκιμασίες, προκρίνοντας σχήματα τρίτης τάξης, που ανήκουν στην οικογένεια σχημάτων TVD και έχουν την ιδιότητα να μεταπίπτουν σε σχήματα μικρότερης τάξης σε περιοχές με έντονη ασυνέχεια, έτσι ώστε να διατηρούν την μονοτονία τους κάτι που δεν υφίσταται με σχήματα, QUICK, SOU και CD. Παράλληλα εξετάζεται ο συνδυασμός σχημάτων χωρικής και χρονικής διακριτοποίησης για την εξαγωγή χρήσιμων συμπερασμάτων με σκοπό την επιλογή του καταλληλότερου συνδυασμού για να περιοριστεί στο ελάχιστο δυνατό το αριθμητικό σφάλμα. Στην συνέχεια, παρουσιάζονται συγκεκριμένες τροποποιήσεις των αριθμητικών εργαλείων που οφείλονται αποκλειστικά στην ΜΥΔ. Μία από αυτές είναι η μέθοδος παρεμβολής για τον υπολογισμό της δύναμης Lorentz. Η ύπαρξη αναλυτικής λύσης, η οποία συγχρόνως είναι ανεξάρτητη από το εξωτερικά επιβαλλόμενο μαγνητικό πεδίο, βοηθά στο να αποδειχθεί ότι η δύναμη Lorentz δεν θα πρέπει να αντιμετωπίζεται ως ένας απλός όρος πηγής των εξισώσεων ορμής, καθώς η κοινή αντιμετώπιση οδηγεί σε αριθμητικά σφάλματα, το μέγεθος των οποίων είναι ανάλογο του τετραγώνου του αριθμού Ha. Παρόλα αυτά όμως η μέθοδος με την καλύτερη ποιότητα αποτελεσμάτων είναι και η δαπανηρότερη σε υπολογιστικό κόστος. Ένα άλλο ενδιαφέρον σημείο που σχετίζεται με την εξίσωση του ηλεκτρικού δυναμικού, βρίσκει εφαρμογή στις προσομοιώσεις τριών διαστάσεων και ιδιαίτερα σε περιπτώσεις με μονωμένα ή πεπερασμένης αγωγιμότητας τοιχώματα. Για την επίλυση των οποίων, απαιτείται η επιπλέον επίλυση της εξίσωσης Poisson για το δυναμικό σε κάθε συνοριακή επιφάνεια, προκειμένου να ικανοποιηθεί η διατήρηση της ρεύματος εντός του υπολογιστικού πεδίου. Οι προσομοιώσεις ΜΥΔ ροής φυσικής συναγωγής αγώγιμού ρευστού αποτελούν ένα προνομιακό χώρο για την πιστοποίηση υπολογιστικών εργαλείων, καθώς υπάρχει πλήθος εργασιών που αποτελούν βάση σύγκρισης, που σχετίζονται με την ενσωμάτωση των ανωτέρω τεχνικών, για να επιτευχθεί αξιόπιστη εκκίνηση για την διερεύνηση άλλων προβλημάτων υπολογιστικής ΜΥΔ. Τα προβλήματα αυτά αφορούν ΜΥΔ ροές σε περιστρεφόμενα εγκλείσματα και η έρευνα εστιάζει στην φυσική αυτών των ροών υπό την επίδραση ομογενούς μαγνητικού πεδίου. Μελετώνται περιπτώσεις κυλινδρικού και σφαιρικού κελύφους, όπου στο εσωτερικό τους υφίσταται αγώγιμο στερεό σώμα που περιστρέφεται με σταθερή γωνιακή ταχύτητα. Στο επίσης αγώγιμο ρευστό αναπτύσσεται παράλληλη δέσμη αντίθετης φοράς από αυτή του περιστρεφόμενου σώματος, η οποία γίνεται λεπτότερη και ταχύτερη με την αύξηση του μαγνητικού πεδίου Τα δύο μεγέθη μεταβάλλονται ασυμπτωτικά με την αύξηση του αριθμού Hartmann. Οι προσομοιώσεις περιορίζονται σε περιοχές αριθμών Ha από 0 μέχρι 600, δίνοντας ιδιαίτερη προσοχή στους μικρούς αριθμούς Ha, καθώς σε αυτές τις περιοχές οι ιξώδεις όροι είναι ισχυροί και υπολογίσιμοι, κάτι που δεν μπορεί να ληφθεί υπόψη στις θεωρητικές προσεγγίσεις. Οι υπολογισμοί σε αυτή την περιοχή επιτρέπουν την εξαγωγή συμπερασμάτων για μικρούς αριθμούς Ha, έτσι ώστε να συμπληρωθεί η εικόνα, που αδυνατούν να παράγουν οι αναλυτικές-ασυμπτωτικές λύσεις στην συγκεκριμένη περιοχή. Στο επόμενο βήμα, εξετάζεται σε κυλινδρική γεωμετρία ποιοτικά η αγωγιμότητα των τοιχωμάτων και οι επιδράσεις αυτών στο πεδίο ροής του αγώγιμου ρευστού. Η διαφορετική αγωγιμότητα των τοιχωμάτων επιφέρει σημαντικές αλλαγές στην συμπεριφορά του ρευστού σε σχέση με την υδροδυναμική περίπτωση. Η κατάλληλη επιλογή των τοιχωμάτων μπορεί να ελέγξει την ροή και την κλίση των ταχυτήτων στην επιφάνεια του τοιχώματος. Επίσης, παρατηρείται ότι η χρήση ηλεκτρικά αγώγιμών τοιχωμάτων οδηγεί στην εμφάνιση υπερ-περιστρεφόμενων δεσμών (super rotational jets) που αποτελεί αφετηρία τεχνολογικών εφαρμογών. Τέλος, στις ενδιάμεσες περιπτώσεις, όπου η ηλεκτρική αγωγιμότητα των τοιχωμάτων είναι μεταξύ των δύο ακραίων περιπτώσεων (πλήρως αγώγιμα ή μονωμένα), η συμπεριφορά του ρευστού δεν είναι αποτέλεσμα γραμμικών υπερθέσεων των δύο άκρων. Αντιθέτως, από τις προσομοιώσεις προκύπτει ότι ακόμη και για πολύ μικρό λόγο αγωγιμότητας, η ροή δεν μπορεί να αντιμετωπιστεί ως ροή με μονωμένα τοιχώματα. Συγχρόνως, ο παραπάνω λόγος αποτελεί βασική παράμετρος για την εξέταση της διάβρωσης των μεταλλικών τοιχωμάτων του εγκλείσματος, με την εφαρμογή απλοποιημένου μοντέλου μεταφοράς συγκεντρώσεων. Με βάση το οποίο παρατηρείται ότι ο ρυθμός διάβρωσης, που σχετίζεται με τον αριθμό Sherwood, μεταβάλλεται μη γραμμικά με το μέτρο αγωγιμότητας.