Προσδιορισμός των υδροδυναμικών διεργασιών της παράκτιας ζώνης με τρισδιάστατη (3D) εφαρμογή του αριθμητικού μοντέλου OpenFOAM και εκτίμηση της αντίστοιχης στερεομεταφοράς με υφιστάμενο αριθμητικό σχήμα

Abstract

Sea waves, propagating large amounts of energy and secondarily mass, from the open sea to the coast, are a very important factor that affects coast stability, the balance of the ecosystem in the coastal zone, the interaction with coastal structures, coastal pollution etc. The appropriate use and exploitation of coastal zones is interconnected with the requirements of detailed planning concerning both technological and manufacturing activities that take place in the coastal environment. This becomes more important taking into account that the coastal environment is an extremely significant factor both in economic and social terms for the overall development of a region. This study focuses on the numerical investigation of hydrodynamic processes and sediment transport rates estimations of coastal zone. Hydrodynamic prediction of significant coastal zone physical wave phenomena such as wave breaking and wave-induced currents is achieved using interfoam multiphase solver, part of OpenFOAM®, a state of the art open-source Computational Fluid Dynamics – CFD, numerical tool. The three dimensional Reynolds – time Averaged Navier - Stokes equations, (RANS) are solved with the contribution of Volume of Fluid method (VOF) for free surface tracking. The numerical scheme is suitably framed by appropriate initial and boundary conditions for each field variable, water-air fraction alpha.water , pressure P_rgh , velocity U , turbulence kinetic Energy (TKE) k, turbulence viscosity nut and dissipation rate of TKE omega. At the inlet and outlet boundaries of the numerical domain, wave generation and absorption conditions were used respectively, as relaxation zones techniques, implemented using waves2Foam, (Jacobsen et al., 2012). The turbulence closure problem is handled using a buoyancy modified k - ω SST model. Specifically, in this work, turbulence modeling is achieved using the incompressible k- ω SST model as implemented by Menter, (1994) with an additional buoyancy term as implemented by Devolder et al., (2018). The buoyancy modified k – ω SST turbulence model differs from the original k – ω SST model as provided in OpenFOAM® , as the density is explicitly included in the turbulence transport equations and the buoyancy term is added to the turbulent kinetic energy (TKE) equation. The buoyancy term according to Devolder et al. (2018), contributes to the suppression of turbulence level at free water surface especially in zones where the governing direction of the density gradient is vertical, i.e. the zone near the free surface where non – breaking waves are propagating and consequently the turbulent viscosity vt, tends to zero. As a result, in case of non-breaking waves the model switches to a laminar regime near the free surface, preventing excessive wave damping. On the contrary, a fully turbulent solution is obtained in the surf zone where the density gradient consists of an important horizontal component. At the breaking point, this condition is obtained when shoaling waves are reaching their limiting wave height. Appropriate interfoam’s setting were used to export successful numerical results as interfoam solver comes with default parameters. The time and space discretization, the iterative solvers options, the solution tolerances and algorithm settings are some of the available numerical parameterization options. In this work it’s worth to refer in this point that in all simulations we use adaptive time stepping based on a maximum allowed CFL number of 0.15 as indicated by Larsen et al., (2019). The pressure-velocity calculation procedure for the Navier-Stokes equations is achieved using the Pressure-Implicit Splitting Operators, (PISO) algorithm. The spatial discretization is performed taking into account at least 10 grid points per wave height and maintaining the ratio 2:1 and 3:1 (depending of each simulation’s characteristics), for the horizontal to vertical dimension of each computational cell.To investigate the accuracy of hydrodynamic numerical results an extensive comparison with benchmark experimental cases was carried out, either in two-dimensional or three-dimensional setup, including the presence of relevant wave phenomena. Specifically, The generation of higher harmonics in waves propagating over submerged obstacles was investigated simulating Beji and Battjes, (1993) work. Results indicate good agreement of wave height distribution, exporting simultaneously, principal characteristics for the spatial discretization of the computational domain. The importance of wave induced currents in surf zone and their impact to sediment transport predictions demands the inspection of such hydrodynamic characteristics of wave flow using the numerical scheme. Initially, a two dimensional (2D) case based on Ting and Kirby, (1994) experimental work was simulated using the combined numerical scheme (interfoam solver and buoyancy modified turbulence model for the prediction of turbulent characteristics of the flow) either for spilling and plunging breakers. The numerical results of the wave flume, concerning wave height distributions, mean water level predictions and variations of time-mean horizontal velocity to depth, indicates good agreement against experimental. Subsequently, the applicability of the numerical scheme was investigated for three dimensional (3D) applications. Firstly, wave diffraction around a detached breakwater submitted to regular wave forcing was examined using a hypothetic numerical wave basin with a breakwater approximately in the middle of the tank. The resulted numerical solution was compared against analytical data of Penny and Price, (1944) work, found in “Wiegel” diffraction table. The comparison achieved by measurements of wave height in the lee of breakwater. The diffraction coefficient, k_d=Η⁄Η_0 as a result of the analytical solution was compared with numerical measurements. The same procedure was followed for the numerical against experimental measurements comparison of wave height distribution, simulating part of Pos and Kilner, (1987) work. To ensure the ability of the numerical scheme to simulate wave induced eddy current pattern behind a detached breakwater, an experimental layout based on Mory and Hamm, (1997) experimental work was investigated concerning wave height and currents measurements in the lee of a detached breakwater under the influence of regular waves. Good consistency between numerical and experimental results is obtained, indicating that the numerical scheme is capable to simulate complex wave induced hydrodynamic processes under the appropriate parameterization. For a complete evaluation and verification of longshore current and relative sediment transport rates, part of the experimental work of Wang et al., (2002) concerning large-scale laboratory measurements of longshore sediment transport under spilling and plunging breakers was simulated. Hydrodynamic results are used to predict total longshore sediment transport rates under the influence of irregular waves. The estimation of total longshore sediment transport rates is achieved using a widespread and widely accepted and used in practice numerical scheme for suspend load and sheet flow transport rate estimation. According to Roelvink and Stive, (1989); Karambas and Koutitas, (2002), time-averaged approach for estimating the suspended load induced by wave breaking can be assumed. In addition, for sheet flow transport rate estimation, Dibajnia and Watanabe, (1992) formula was used. Thus, the total submerged weight transport rate (including additional stirring of sediment by the surface breaking-induced turbulence that penetrates toward the bottom) can be estimated using the equation 5.18. Finally, to ensure the applicability of the numerical scheme as a tool in realistic dimensions coastal engineering application concerning the prediction of the hydrodynamic characteristic of the flow, a numerical wave basin of 148-m cross-shore, 340-m longshore dimensions and maximum water depth of 5-m, was simulated under the influence of irregular waves. Exported results concerning cross-shore distributions of longshore current and turbulence kinetic energy were used to predict total longshore sediment transport rates. The Kamphuis (1991) empirical formula was used to validate exported numerical results. Concisely, the present work demonstrates the effectiveness of the open source (CFD) computational package, OpenFOAM®, by assessing the hydrodynamic processes of coastal zone. The derived hydrodynamic results concerning wave height, turbulent kinetic energy and longshore current are used to predict longshore sediment transport rates. The combined numerical scheme aims to demonstrate an open source, scientifically validated powerful numerical tool, including all stages of a numerical simulation (pre-processing, processing, post processing) for coastal zone’s hydrodynamic predictions and coastal technical structures design. Concurrently, this work establish the base for a full three dimensional hydro-morphodynamic numerical scheme development for relative coastal zone applications.

Οι διαδιδόμενοι από τα ανοιχτά προς την ακτή θαλάσσιοι κυματισμοί, μεταφέρουν πρωτίστως μεγάλα ποσά ενέργειας και δευτερευόντως μάζας, επιδρώντας άμεσα στην ισορροπία του παράκτιου περιβάλλοντος. Η διασφάλιση της εύρυθμης λειτουργιάς της παράκτιας ζώνης συνάδει με τις έννοιες της κοινωνικής και οικονομικής ευημερίας και ανάπτυξης, καθώς πλήθος ανθρωπογενών δραστηριοτήτων αναπτύσσονται στη περιοχή αυτή. Ένα από τα βασικά αίτια υποβάθμισης του παρακτίου περιβάλλοντος είναι αυτό της διάβρωσης των ακτών, διεργασία η οποία προέρχεται είτε από επιβλαβείς ανθρωπογενείς παρεμβάσεις (άστοχα τεχνητά έργα, κ.α.) είτε από φυσικές διεργασίες όπως η εμφάνιση της κλιματικής αλλαγής. Στην κατεύθυνση αυτή, η παρούσα διατριβή επιχειρεί να δώσει μια ολοκληρωμένη απάντηση στο πρόβλημα της διάβρωσης του παρακτίου μετώπου η οποία βασίζεται στη χρήση προηγμένων μαθηματικών ομοιωμάτων παράκτιας υδροδυναμικής και προσδιορισμού της στερεομεταφοράς κατά μήκος των ακτών. Για τον προσδιορισμό των υδροδυναμικών χαρακτηριστικών της κυματικής ροής χρησιμοποιείται ο επιλυτής πολλαπλών φάσεων interfoam οποίος περιλαμβάνεται στο πακέτο υπολογιστικής ρευστοδυναμικής (CFD) ανοιχτού κώδικα OpenFOAM® , σε συνδυασμό με τη κυματική βιβλιοθήκη waves2Foam, (Jacobsen et al. 2012). H επίλυση των τρισδιάστατων μέσων χρονικά εξισώσεων ροής Navier – Stokes, (RANS) επιτυγχάνεται με τη συνεισφορά της μεθόδου Volume Of Fluid, (VOF), για την ανίχνευση της διεπιφάνειας των δύο ρευστών, αέρας – νερό, στα οπoία διαδίδεται η κυματική διαταραχή. Το αριθμητικό σχήμα εμπλουτίζεται με κατάλληλες αρχικές και συνοριακές συνθήκες (οι οποίες παρέχονται μέσω του OpenFOAM®) για κάθε μεταβλητή του πεδίου ροής, όπως το κλάσμα αέρα νερού alpha.water σε κάθε κελί του υπολογιστικού χωρίου, η πίεση P_rgh (Ολική πίεση μείον υδροστατική), η ταχύτητα των σημείων του πεδίου ροής, U , η τυρβώδης κινητική ενέργεια ΤΚΕ, k , το τυρβώδες ιξώδες, nut, καθώς και ο ρυθμός σκέδασης της τυρβώδους κινητικής ενέργειας omega. Στην είσοδο και στην έξοδο του υπολογιστικού χωριού εφαρμόζονται ζώνες παραγωγής και απορρόφησης κυματισμού οι οποίες εφαρμόζονται με χρήση της κυματικής βιβλιοθήκης waves2Foam. Η προσομοίωση των τυρβωδών χαρακτηριστικών της κυματικής ροής και συνεπώς η διαχείριση του προβλήματος «κλεισίματος της τύρβης» στην παρούσα εργασία επιτυγχάνεται με εφαρμογή του μοντέλου τύρβης k – ω SST όπως αυτό αναπτύχτηκε από τον Menter (1994). Το συγκεκριμένο μοντέλο χρησιμοποιείται στην παρούσα περίπτωση αφενός με την προσθήκη του όρου της πυκνότητας στις δυο εξισώσεις του, ούτως ώστε να εναρμονίζεται με τη λειτουργία του επιλυτή δύο φάσεων interfoam και τη Μέθοδο VOF. Αφετέρου, με μια προσθήκη η οποία προέρχεται από την εργασία των Devolder et al. (2018) και στοχεύει στην απαλοιφή της τυρβώδους κινητικής ενέργειας που εμφανίζεται στη διεπιφάνεια νερού αέρα σε περιοχές στρωτής ροής, όπου συμφώνα με τους Larsen and Fuhrman, (2018) , δεν δικαιολογείται η εμφάνιση της. Αυτό επιτυγχάνεται με την εισαγωγή ενός επιπρόσθετου όρου Gb (buoyancy term) στην εξίσωση της τυρβώδους κινητικής ενέργειας k, ο οποίος ουσιαστικά αφαιρεί το παραγόμενο ποσό τυρβώδους κινητικής ενέργειας στις περιοχές στρωτής ροής. H επιλογή παραμέτρων του επιλυτή interfoam ούτως ώστε να είναι σε θέση να αποτυπώσει με ακρίβεια τα χαρακτηριστικά της κυματικής ροής για τα οποία εφαρμόζεται, επιτυγχάνεται στην παρούσα εργασία μέσα από συμπεράσματα που προέκυψαν κατά την εφαρμογή, όπως ο βέλτιστος αριθμός σημείων ανά ύψος κύματος, καθώς και με υποδείξεις από τη συναφή βιβλιογραφία, π.χ. Larsen et al., (2019). Για την ικανοποίηση του κριτηρίου Courant χρησιμοποιήθηκε η υπόδειξη των Larsen et al., (2019), οι οποίοι χρησιμοποιώντας τον επιλυτή interfoam για τον προσδιορισμό της διάδοσης μοναχικού κύματος σε περιοχή σταθερού βάθος νερού, χρησιμοποίησαν το κριτήριο Courant, CFL ≤ 0.15, καταγράφοντας το φαινόμενο με ακρίβεια. Επιπρόσθετα οι εξισώσεις ροής επιλύονται με χρήση του αλγορίθμου Pressure-Implicit with Splitting of Operators, (PISO), σε όλες τις προσομοιώσεις, ενώ τέλος από παρατηρήσεις που προέκυψαν κατά την διεξαγωγή των επιμέρους προσομοιώσεων, η ακρίβεια στην καταγραφή του ύψους κύματος φαίνεται να απαιτεί 10 σημεία υπολογιστικού πλέγματος ανά ύψος κύματος κατά την αρχική διαταραχή, ορίζοντας έτσι την κατακόρυφη διακριτοποίηση του υπολογιστικού πλέγματος. Οι οριζόντιες διαστάσεις δε, διακριτοποιήθηκαν συγκριτικά με την κατακόρυφη διάσταση χρησιμοποιώντας λόγο 2:1 ή 3:1 (οριζόντια προς κατακόρυφη διάσταση). Κατά την εφαρμογή του αριθμητικού σχήματος, αρχικά διερευνήθηκε η αποτελεσματικότητα του στην αποτύπωση βασικών κυματικών διεργασιών μέσα από τη προσομοίωση σημαντικών πειραματικών εργασιών στις οποίες συναντώνται τα προς μελέτη κυματικά φαινόμενα. Τα εξαγόμενα αριθμητικά αποτελέσματα συγκρίνονται σε κάθε περίπτωση με τα αντίστοιχα πειραματικά. Η διαμόρφωση του ύψους κύματος σε κανάλι μεταβαλλόμενης βυθομετρίας με βυθισμένο ύφαλο μελετήθηκε μέσα από τη προσομοίωση του πειράματος των Beji and Battjes, (1993). Η παραγωγή αρμονικών κυματομορφών ανώτερης τάξης φάνηκε να καταγράφεται με ακρίβεια από το αριθμητικό σχήμα ενώ παράλληλα, μέσα από τη συγκεκριμένη εφαρμογή προέκυψαν βασικά συμπεράσματα σχετικά με τη χωρική διακριτοποίηση του αριθμητικού πλέγματος. Βασικό παράγοντα διαμόρφωσης του παράκτιου μετώπου μέσω της στερεομεταφοράς αποτελεί η κυκλοφορία κυματογενούς προέλευσης. Σε πρώτη φάση διερευνήθηκε η ικανότητα του αριθμητικού ομοιώματος στην αποτύπωση υδροδυναμικών διεργασιών και συγκεκριμένα στη ορθή καταγραφή του υποβρύχιου ρεύματος επαναφοράς, undertow. Τα αριθμητικά αποτελέσματα σχετικά την κατακόρυφη κατανομή των ταχυτήτων στο πεδίο ροής δυσδιάστατου καναλιού μεταβαλλόμενης βυθομετρίας συγκριθήκαν με τα αντίστοιχα πειραματικά της εργασίας των Ting and Kirby, (1994).Σε επόμενο στάδιο, με τα συμπεράσματα που προέκυψαν από τις δισδιάστατες – χωρικά προσομοιώσεις συντελείται η τρισδιάστατη εφαρμογή του αριθμητικού ομοιώματος εισάγοντας μια ακόμη οριζόντια διάσταση, γεγονός που αυξάνει κατακόρυφα της απαίτηση υπολογιστικών πόρων, αλλά ταυτόχρονα και τη δυνατότητα του αριθμητικού ομοιώματος στη καταγραφή σύνθετων κυματικών και υδροδυναμικών διεργασιών. Συγκεκριμένα, με σκοπό τη μελέτη και καταγραφή του φαινομένου της περίθλασης διαδιδόμενων επιφανειακών κυματισμών από το άκρο ημιάπειρου κυματοθραύστη δημιουργήθηκε αντίστοιχη υποθετική δεξαμενή στη οποία λαμβάνει χώρα το σχετικό φαινόμενο. Τα αριθμητικά αποτελέσματα του συντελεστή περίθλασης k_d=Η⁄Η_0 , όπου Η το ύψος κύματος σε συγκεκριμένες θέσεις στη «σκιά» του κυματοθραύστη και H0 το ύψος κύματος στην κορυφή (είσοδο) του κυματοθραύστη, συγκρίνονται με την αναλυτική λύση των Penny and Price, (1944). Παράλληλα μελετήθηκε η σύγκριση αντίστοιχων αριθμητικών αποτελεσμάτων με αντίστοιχα πειραματικά μέσα από την προσομοίωση μέρους της πειραματικής εργασίας των Pos and Kilner, (1987). Η ικανότητα του αριθμητικού σχήματος στην αποτύπωση σύνθετων υδροδυναμικών διεργασιών όπως ρεύματα υπό μορφή στροβίλων, διερευνήθηκε μέσα από τη προσομοίωση της πειραματικής εργασίας των Mory and Hamm (1997) όπου κατά την κάθετη πρόσπτωση επιφανειακών αρμονικών κυματισμών σε κυματοθραύστη δημιουργείται ισχυρός στρόβιλος στη σκιά αυτού, με σχεδόν αδιατάρακτο κέντρο περιστροφής και υψηλές ταχύτητες στην περιφέρεια του, υποδεικνύοντας καλή συμφωνία με τις αντίστοιχες πειραματικές μετρήσεις. Σε επόμενο στάδιο, με γνώμονα τη βασική συνεισφορά του παράλληλου στη ακτογραμμή ρεύματος στη στερεοπαροχή της παράκτιας ζώνης, προσομοιώθηκε μέρος της πειραματικής εργασίας των Wang et al., (2002), κατά την οποία καταγράφεται η εγκάρσια κατανομή των ταχυτήτων του διαμήκους ρεύματος το οποίο δημιουργείται από την πρόσπτωση τυχαίων κυματισμών υπό γωνία σε ακτή ήπιας κλίσης. Τα αριθμητικά αποτελέσματα αναδεικνύουν την ύπαρξη μέγιστων τιμών ρεύματος ακόμη και λίγο πριν τη ζώνη αναρρίχησης του κυματισμού στην ακτή. Κατά την εκτίμηση της συνολικής στερεοπαροχής, για το φορτίο σε αιώρηση που δημιουργείται από τη θραύση του κύματος, χρησιμοποιήθηκαν μέσες χρονικά προσεγγίσεις εφαρμόζοντας τη μέθοδο Baillard (1991). Η συνολική δε στερεοπαροχή προσδιορίστηκε λαμβάνοντας υπόψη το φορτίο πυθμένα σύμφωνα με τους Dibajnia and Watanabe, (1992) και το φορτίο σε αιώρηση σύμφωνα με τους Roelvink and Stive (1989); Karambas and Koutitas (2002). Τα εξαγόμενα αριθμητικά αποτελέσματα έρχονται σε συμφωνία με τα αντίστοιχα πειραματικά τόσο για την περίπτωση θραύσης εκχείλισης όσο και για την περίπτωση θραύσης εκτίναξης. Σε τελικό στάδιο διερευνήθηκε η εφαρμοσημότητα του αριθμητικού σε υποθετική ακτή ρεαλιστικών διαστάσεων 340m μήκος, 148m πλάτος και μέγιστο βάθος νερού 5m (scale up), ήπιας κλίσης. Για διάφορες τιμές ύψους και περιόδου κύματος, μελετήθηκε η ανάπτυξη του παράλληλου στην ακτή ρεύματος υπό διαφορετικές γωνίες πρόσπτωσης φάσματος κυματισμού τύπου JONSWAP, στην ακτή. Η εκτιμώμενη στερεομεταφορά που προέκυψε επαληθεύτηκε με χρήση της εμπειρική φόρμουλας του Kamphuis (1991). Συνοπτικά μπορεί να ειπωθεί ότι κατά την παρούσα διατριβή μελετήθηκε η αποτελεσματικότητα του αριθμητικού ομοιώματος που βασίζεται στον επιλυτή πολλαπλών φάσεων interfoam στην αποτίμηση των κυματικών και υδροδυναμικών διεργασιών της παράκτιας ζώνης. Τα εξαγόμενα αποτελέσματα χρησιμοποιούνται για τον προσδιορισμό της στερεοπαροχής με απώτερο σκοπό την ανάδειξη ενός σύγχρονου ισχυρού αριθμητικού μοντέλου για την εφαρμογή σε συναφή προβλήματα ακτομηχανικής. Ταυτόχρονα, δημιουργείται η βάση, για τη δόμηση ενός ολοκληρωμένου υδρο-μορφοδυναμικού μοντέλου στη παράκτια ζώνη, για το σχεδιασμό έργων προστασίας ακτών και την εκτίμηση των επιπτώσεων των παράκτιων τεχνικών έργων στη μορφολογία της ακτής.