Topology and shape optimization in fluid mechanics & conjugate heat transfer using continuous adjoint with consistent discretization

Abstract

This PhD thesis develops methods to perform adjoint-based Topology Optimization (TopO) in Computational Fluid Dynamics (CFD) problems with (or without) Conjugate Heat Transfer (CHT). All methods and tools are developed within the open–source CFD code of OpenFOAM, extending the adjointOptimisation library programmed and made publicly available by the Parallel CFD & Optimization Unit of the NTUA. A key contribution of this thesis is to alleviate some of the well–known weaknesses of the widely–used density–based TopO (denTopO) approach, related to the accuracy of the corresponding analysis code. denTopO is performed on a static CFD grid by introducing a single design variable per grid cell which stands for an impermeability and indicates whether each cell is occupied by fluid or solid material. Ideally, binary designs, where all the design variables take on 0/1 values, are expected. In view of this and to avoid well–known checkerboard pathologies, the filtering of the design variables based on a Helmholtz PDE is performed. This yields a smoother field that is later projected using a smooth/differentiable Heaviside function to become sharper. The so-computed field is introduced to the flow equations through the Brinkman penalization terms, standing for blockage forces inside the solid parts of the domain. The flow is modeled by the Reynolds-Averaged Navier–Stokes (RANS) equations coupled with the Spalart–Allmaras turbulence model. CHT problems require solving the temperature equation over the fluid and solid domains. Due to the one–way coupling between the mean flow and the temperature equation, the latter is solved as a post–processing step by making use of proper interpolation schemes for the thermophysical properties based on the values of the design variables. This thesis relies upon the continuous adjoint method which computes the gradients of the objective and constraint functions with respect to (w.r.t.) the design variables. This computation comes at a cost that does not depend on the number of design variables which, in TopO, is equal the number of grid cells.The ability of the denTopO method to compute innovative designs and unexpected topologies comes with lower accuracy of the CFD/CHT solution, as the denTopO solver is unable to impose exact conditions on the Fluid Solid Interfaces (FSIs). Instead, it models solid areas as porous regions of high impermeability, without being able to avoid spurious/unwanted flow leakage. Extracting the FSI from the TopO solution is a task that, depending on the way it is performed and the formulation of TopO (e.g. density–based, level–set etc.), might lead to slightly different optimized solutions. For this reason, and due to the aforementioned inaccuracies of denTopO, performances may change if the optimized solutions from denTopO are re–evaluated on body–fitted grids. As a first remedy, this thesis proposes to compute the FSI upon completion of denTopO and, then, generate body–fitted grids in order to proceed with Shape Optimization (ShpO) as a means to further refine the denTopO solutions; however, this increases the computational cost. To eliminate the lack of accuracy of the denTopO solver in the near wall regions and avoid continuing with ShpO, this thesis proposes and implements the cut–cell method which computes the intersections of the FSI with the cells of the background CFD grid in each cycle of TopO. This leads to the cut–cell TopO method which is still based on the impermeability field of denTopO. However, in contrast to denTopO, this is used exclusively to compute the FSI. The solution of the flow and CHT equations is, thus, performed, on the cut–cell grid, allowing the imposition of exact boundary conditions on the FSI. Moreover, cut–cell TopO applies h-refinement close to the computed FSIs. ShpO is no more needed, and accurate performance values as well as a clear FSI become available. To update the design variables, the derivatives of the objective and constraints w.r.t. the node positions on the FSIs are computed first, and are then transformed into derivatives of the same variables w.r.t. cell impermeabilities by applying the chain rule of differentiation. The superiority of the cut–cell TopO method w.r.t. the ‘‘standard’’ denTopO method is assessed by comparing the results in different test cases. The next part of this dissertation pertains to the TopO of two-fluid heat exchangers, performed by extending the previously developed denTopO and cut–cell TopO tools. In such cases, it is important to ensure the non-mixing of the two fluids. The developed framework relies on a single field of design variables to parameterize regions occupied by solid material or either of the two working fluids. Compared to mono–fluid TopO, the bi–fluid TopO method incorporates the filtering and projection operators in two steps. The second step controls the user-specified minimum solid thickness separating the two fluids, and yields two distinct fluid indicator fields. In denTopO, these are used to impose the Brinkman terms in the flow equations, solved separately for each fluid. In cut–cell TopO, the two fluid indicators are used to compute the FSI and generate cut–cells, without affecting directly the governing equations. The developed tools for TopO, either based on the ‘‘standard’’ density approach or the cut–cell method, are assessed in problems of academic and of practical interest. 2D internal flow problems targeting min. total pressure drop and/or max. heat exchange are tackled. The application of mono–fluid TopO for the minimization of total pressure losses inside an HVAC duct, bringing cool air from the front console to the rear part of a passenger car, is also worked out. For bi–fluid TopO, both the density-based and the cut-cell method are assessed in the optimization of a compact heat exchanger with many inlets and outlets, for the two working fluids.As mentioned before, both ShpO and TopO are supported by continuous adjoint. Selecting proper discretization schemes for the solution of the continuous adjoint problem affects the accuracy of the computed derivatives. On the other hand, discrete adjoint solvers provide consistent gradients, as the discretized primal residuals are differentiated. This PhD thesis develops the ‘‘Think Discrete-Do Continuous’’ (TDDC) adjoint for the incompressible RANS equations solved in OpenFOAM using the SIMPLE algorithm for the pressure–velocity coupling. In the TDDC adjoint, the development starts by hand–differentiating the discretized residuals of the primal solver, as in discrete adjoint. The terms emerging from this development are appropriately re–arranged to reverse–engineer consistent fluxes for the discretization of all terms in the continuous adjoint PDEs. The new adjoint solver combines the best of discrete and continuous adjoint, as it gives exact/consistent gradients of the objective function which are in perfect agreement with those computed by Finite Differences (FDs) and maintains a low memory footprint. This dissertation provides two ways of solving the set of adjoint PDEs. The first is based on a fixed–point iterative scheme derived by first rewriting the segregated SIMPLE algorithm in the form of a fully–coupled left–preconditioned Fixed Point Iteration (FPI). The transpose of the derived preconditioning matrix is then used for the adjoint FPI. Even in fully coupled form, the preconditioner finally acts on the adjoint residuals in a segregated manner. This leads to the adjoint to the SIMPLE algorithm that performs the exact same operations as the primal algorithm but in reverse order. This strategy is duality–preserving as the adjoint solver inherits the convergence characteristics of the primal FPI. However, in real–world applications the primal solver is likely not to converge tightly/monotonically, but enter Limit Cycle Oscillations (LCO). In such cases, solving the steady–state adjoint problem by making a linearization around the last snapshot of the flow solution is likely to diverge since the adjoint FPI (derived directly from the primal FPI) might not be contractive, i.e. might possess eigenvalues of magnitude greater than 1. To address this issue, the solution of the adjoint problem in a fully coupled manner by making use of Krylov subspace linear solvers is proposed. The full adjoint system is solved using the Flexible variant of the Generalized Minimal Residual method with Deflated Restarting (FGMRES-DR) between restarts to avoid residual stagnation. The linear solver is implemented in OpenFOAM and makes use of the developed adjoint SIMPLE algorithm as preconditioner. The coupled Krylov adjoint solver is robust enough to converge even in problems where the primal solver enters LCO, with enough, though, Krylov basis vectors. Occasionally, the coupled adjoint solver may significantly reduce the cost of solving the adjoint problem, at the expense of a larger memory footprint.

Η Διδακτορική αυτή Διατριβή πραγματεύεται την ανάπτυξη, πιστοποίηση και εφαρμογή μεθόδων και λογισμικού για προβλήματα βελτιστοποίησης μορφής και κυρίως τοπολογίας τόσο στη ρευστοδυναμική όσο και σε προβλήματα Συζευγμένης Μεταφοράς Θερμότητας (ΣΜΘ). Οι μέθοδοι αναπτύσσονται στο ανοιχτό λογισμικό OpenFOAM επεκτείνοντας τη βιβλιοθήκη adjointOptimisation που έχει αναπτυχθεί, και είναι προσβάσιμη από όλους, από τη Μονάδα Παράλληλης Υπολογιστικής Ρευστοδυναμικής & Βελτιστοποίησης (ΜΠΥΡ&Β) του ΕΜΠ. Σημαντική συνεισφορά αυτής της Διατριβής είναι η επέκταση λογισμικού Βελτιστοποίησης Τοπολογίας (ΒΤοπ) που αναπτύχθηκε σε προηγούμενες διδακτορικές διατριβές της ΜΠΥΡ&Β για προβλήματα Ρευστοδυναμικής σε προβλήματα ΣΜΘ τόσο με ένα όσο και με δύο ρευστά εμποδίζοντας παράλληλα την ανάμιξή τους. Γίνεται χρήση της γνωστής στη βιβλιογραφία μεθόδου της τεχνητής αδιαπερατότητας η οποία εισάγει μια μεταβλητή σχεδιασμού ανά όγκο ελέγχου του υπολογιστικού πλέγματος η οποία και υποδηλώνει την παρουσία ρευστού ή στερεού εντός του όγκου ελέγχου. Για την εξάλειψη ανεπιθύμητων “μοτίβων σκακιέρας” κατά τη ΒΤοπ πραγματοποιείται ομαλοποίηση του πεδίου αδιαπερατότητας μέσω της επίλυσης μιας μερικής διαφορικής εξίσωσης (ΜΔΕ) τύπου Helmholtz. Έπειτα, γίνεται χρήσης μιας ομαλής/διαφορίσιμης συνάρτησης Heaviside ούτως ώστε το τελικό πεδίο αδιαπερατότητας να λαμβάνει τιμές 0/1 που υποδηλώνουν ρευστό/στερεό, αντίστοιχα. Το πεδίο αυτό εισάγεται στις εξισώσεις ασυμπίεστης ροής, που επιλύονται με τον αλγόριθμο SIMPLE, μέσω των όρων ποινής Brinkman ώστε οι ροϊκές μεταβλητές να λαμβάνουν (πρακτικά) μηδενικές τιμές εντός του στερεού χωρίου. Το αναπτυχθέν λογισμικό εφαρμόζεται τόσο σε προβλήματα Ρευστοδυναμικής όσο και ΣΜΘ. Στη Ρευστοδυναμική, όπου η ροή είναι τυρβώδης και χρησιμοποιείται το μοντέλο τύρβης των Spalart-Allmaras, γίνεται βελτιστοποίηση ενός 2Δ τετραγωνικού χωρίου σχεδιασμού με μία είσοδο και δύο εξόδους για το ρευστό και ενός 3Δ αγωγού κλιματισμού επιβατικού αυτοκινήτου. Σε προβλήματα ΣΜΘ με ένα ρευστό, γίνεται βελτιστοποίηση της ψύξης 2Δ πηγών θερμότητας σε διαφορετικές διατάξεις και πραγματοποιούνται παραμετρικές μελέτες για την επίδραση της ομαλοποίησης στη βέλτιστη λύση. Ακόμα, πραγματοποιείται ΒΤοπ ενός 3Δ εναλλάκτη θερμότητας τύπου αντιροής με τα δύο ρεύματα να εισέρχονται/εξέρχονται από τον εναλλάκτη από πολλαπλές εισόδους/εξόδους, επιβάλλοντας ίση κατανομή της παροχής του κρύου ρευστού σε όλες τις εξόδους κατά τη ΒΤοπ. Το κύριο πλεονεκτήματα της ΒΤοπ με τη μέθοδο της τεχνητής αδιαπερατότητας είναι η δυνατότητα να αλλάζει δραστικά την τοπολογία της λύσης καθώς και η δυνατότητα υπολογισμού μη-συμβατικών λύσεων η οποία οφείλεται στο μεγάλο πλήθος μεταβλητών σχεδιασμού. Αυτά συνοδεύονται συνήθως από ελάχιστο κόπο όσον αφορά την ανάπτυξη του σχετικού λογισμικού. Το μειονέκτημα της μεθόδου είναι η μειωμένη ακρίβεια υπολογισμού της ροής καθότι δεν υπολογίζεται η Διεπιφάνεια Ρευστού Στερεού (ΔΡΣ) με αποτέλεσμα η επιβολή των οριακών συνθηκών του τοίχου να γίνεται έμμεσα μέσω των όρων ποινής τύπου Brinkman. Για αυτόν τον λόγο, οι λύσεις της ΒΤοπ επαναξιολογούνται σε οριόδετα πλέγματα, αφότου γίνει υπολογισμός της ΔΡΣ, με την επαναξιολόγηση να οδηγεί συνήθως στον υπολογισμό διαφορετικών ολοκληρωματικών μεγεθών. Σε αυτήν τη Διατριβή, προτείνονται δύο εναλλακτικοί τρόποι για την αντιμετώπιση συναφών ανακριβειών της ΒΤοπ: (α) ο υπολογισμός της ΔΡΣ στο τέλος της ΒΤοπ και η γένεση οριόδετου πλέγματος ώστε να ακολουθήσει ένα επιπλέον στάδιο σχεδιασμού όπου πραγματοποιείται Βελτιστοποίηση Μορφής (ΒΜορφ), (β) η ανάπτυξη μιας εξ ολοκλήρου νέας μεθόδου ΒΤοπ η οποία βασίζεται στη Μέθοδο των Τεμνόμενων Κυψελών (ΒΤοπΤΚ). Η ΒΤοπΤΚ κάνει χρήση του πεδίου τεχνητής αδιαπερατότητας ώστε βάσει αυτού να ορίσει τη ΔΡΣ την οποία και υπολογίζει σε κάθε κύκλο βελτιστοποίησης. Με αυτόν τον τρόπο, γίνεται επίλυση των εξισώσεων ροής σε οριόδετο πλέγμα με επιβολή σωστών οριακών συνθηκών στα όρια ρευστού-στερεού. Ταυτόχρονα, η ΒΤοπΤΚ πραγματοποιεί Προσαρμοστική Πύκνωση του Πλέγματος (ΠΠΠ) που αποσκοπεί όχι στην αύξηση του πλήθους των μεταβλητών σχεδιασμού αλλά στην περαιτέρω βελτίωση της ακρίβειας του ροϊκού επιλύτη. Η σημασία χρήσης ΠΠΠ κατά τη ΒΤοπΤΚ επιβεβαιώνεται με παραμετρικές μελέτες. Τέλος, επιβεβαιώνεται η υπεροχή της νέας μεθόδου αναφορικά με την “κλασική” μέθοδο ΒΤοπ που βασίζεται στην επιβολή των όρων τύπου Brinkman στις εξισώσεις συγκρίνοντας τα αποτελέσματα των δύο μεθόδων σε 2Δ και 3Δ προβλήματα.Η ΒΤοπ, είτε κάνοντας χρήση των όρων Brinkman είτε της μεθόδου των τεμνόμενων κυψελών, υποστηρίζεται από τη συνεχή συζυγή μέθοδο για τον υπολογισμό των παραγώγων ευαισθησίας οι οποίες απαιτούνται για την ανανέωση των μεταβλητών σχεδιασμού. Το κύριο πλεονέκτημα της συνεχούς συζυγούς μεθόδου είναι η φυσική εποπτεία των συζυγών ΜΔΕ και των οριακών συνθηκών τους, και, ειδικότερα, οι χαμηλές απαιτήσεις σε μνήμη του σχετικού λογισμικού. Τα σχήματα διακριτοποίησης που χρησιμοποιούνται για τις συζυγείς εξισώσεις μπορούν να επηρεάσουν σημαντικά την ακρίβεια υπολογισμού των παραγώγων ευαισθησίας, όσο αυτά δεν είναι συμβατά με το διακριτοποιημένο ροϊκό (πρωτεύον) πρόβλημα. Η διακριτή συζυγής μέθοδος εξασφαλίζει αυτή τη συμβατότητα, αλλά συνήθως με αυξημένες απαιτήσεις μνήμης. Η Διατριβή γεφυρώνει το χάσμα ανάμεσα στις δύο συζυγείς μεθόδους αναπτύσσοντας μια παντελώς νέα “Think Discrete-Do Continuous” (TDDC) συζυγή μέθοδο στις ασυμπίεστες εξισώσεις ροής που επιλύονται με την τεχνική διόρθωσης της πίεσης (αλγόριθμος SIMPLE). H TDDC συζυγής μέθοδος εμπνέεται από τη διακριτή συζυγή μέθοδο προκειμένου να αναπτύξει συμβατά σχήματα διακριτοποίησης για της συνεχείς συζυγείς ΜΔΕ. Η συμβατή διακριτοποίηση που επιτυγχάνει η νέα μέθοδος υπερτερεί της “κλασικής” ως προς την ακρίβεια υπολογισμού των παραγώγων ευαισθησίας, συγκρίνοντας με τις παραγώγους που υπολογίζει η μέθοδος των Πεπερασμένων Διαφορών, ανεξαρτήτως της πυκνότητας του πλέγματος που χρησιμοποιείται για τον υπολογισμό αυτών. Αυτό πιστοποιείται σε ένα πρόβλημα εσωτερικής αεροδυναμικής όπου μελετάται ένας 2Δ αγωγός τύπου S. Η βελτιωμένη ακρίβεια υπολογισμού παραγώγων ευαισθησίας συντελεί και στην καλύτερη πορεία ενός αλγορίθμου αιτιοκρατικής βελτιστοποίησης.Τέλος, η Διατριβή προτείνει δύο εναλλακτικούς και αποδοτικούς τρόπους για την επίλυση των συζυγών εξισώσεων. Ο πρώτος βασίζεται στην ανάπτυξη ενός επαναληπτικού σχήματος μέσω του οποίου επιτυγχάνεται δυϊκότητα του συζυγούς επιλύτη με την έννοια ότι ο τελευταίος κληρονομεί τα χαρακτηριστικά και τον ρυθμό σύγκλισης του ροϊκού (πρωτεύοντος) επιλύτη. Η επαναληπτική διαδικασία, επονομαζόμενη “ο συζυγής αλγόριθμος SIMPLE”, αναπτύσσεται αντιστρέφοντας τα βήματα του πρωτεύοντος αλγορίθμου και ενδείκνυται για χρήση σε προβλήματα όπου ο ροϊκός επιλύτης συγκλίνει μονότονα (χωρίς τα υπόλοιπα της ροής να ταλαντώνονται) καθώς εξασφαλίζει τη σύγκλιση των συζυγών εξισώσεων με τον ίδιο ρυθμό. Ωστόσο, σε πολλές εφαρμογές, ο ροϊκός επιλύτης δεν συγκλίνει πλήρως. Αντ’ αυτού τα υπόλοιπα των εξισώσεων ροής τείνουν να ταλαντώνονται. Σε αυτές τις περιπτώσεις, η χρήση του δυϊκού συζυγούς επιλύτη μπορεί να οδηγήσει στην απόκλιση των συζυγών εξισώσεων αν η γραμμικοποίηση πραγματοποιηθεί γύρω από ένα ροϊκό στιγμιότυπο για το οποίο το μητρώο του πρωτεύοντος επαναληπτικού σχήματος έχει ιδιοτιμές μεγαλύτερες της μονάδας. Για τη σταθεροποίηση του συζυγούς προβλήματος σε τέτοιες περιπτώσεις, η Διατριβή αυτή προτείνει την επίλυση των συζυγών εξισώσεων με πλήρως πεπλεγμένο τρόπο κάνοντας χρήση μεθόδων υπόχωρων Krylov. Συγκεκριμένα, γίνεται προγραμματισμός της ευέλικτης εκδοχής του γραμμικού επιλύτη GMRES (FGMRES). Προκειμένου να ενισχυθεί η πορεία σύγκλισης του γραμμικού επιλύτη και να είναι συγκρίσιμη με εκείνη της πλήρους FGMRES, γίνεται υπολογισμός προσεγγιστικών ιδιοτιμών και ιδιοδιανυσμάτων του Ιακωβιανού μητρώου για την επανεκκίνηση της διαδικασίας Arnoldi.