The Extended Finite Element Method for Crack Propagation Problems: Theory and Implementation Details

Μεταπτυχιακός Φοιτητής : Μπακαλάκος Σεραφείμ              
Επιβλέπων Καθηγητής: Παπαδρακάκης M., Καθηγητής
Ημερομηνία : Οκτώβριος 2017

 Η εκτετεμάνη μέθοδος πεπερασμένων στοιχείων (eXtended Finite Element Method – XFEM) αρχικά αναπτύχθηκε από τους Belytsckho, Moes  και Dolbow για την αντιμετώπιση προβλημάτων που αφρούν την διάδοση ρωγμών. Στα προβλήματα αυτά το πεδίο μετατοπίσεων είναι ασυνεχές εγκάρσια στη ρψγμή και το πεδίο τάσεων ιδιάζον στην αιχμή της ρωγμής. Για την επίλυση τους με την παραδοσιακή μέθοδο πεπερασμένων στοιχείων (FEM), απαιτείται η χρήση ενός άριθμητικού πλέγματος που να προσαρμόζεται στη γεωμετρία της ρωγμής, έτσι ώστε κανέναν στοιχείο να μην τέμνεται από αυτή. Κάθε φορά που η ρωγμή διαδίδεται, το πλέγμα αυτό θα πρέπει να επαναδημιουργείται καθιστώντας τη διαδικασία χρονοβόρα για πολύπλοκες γεωμετρίες φορέων. Επιπλέον κάθε φορά θα πρέπει να απεικονίζεται το πεδίο μετατοπίσεων μεταξύ του παλιού και νέου πλέγματος πεπερασμένων στοιχείων, εισάγοντας έτσι σφάλματα στην επίλυση.  

Αντίθετα η XFEM χρησιμοποιεί ένα σταθερό πλέγμα πεπερασμένων στοιχείων, αλλά το πεδίο μετατοπίσων που αναζητεί μπορεί να είναι ασυνεχές. Αυτό επιτυγχάνεται εμπλουτίζοντας τις πολυωνυμικές συναρτήσεις βάσης της παραδοσιακής FEM με ειδικές συναρτήσεις εμπλουτισμού. Για την μοντελοποίηση του άλματος στο πεδίο μετατοπίσεων, χρησιμοποιούμε ως συνάρτηση εμπλουτισμού την βηματική συνάρτηση Heaviside, πολλαπλασιάζοντάς την με τις συμβατικές συναρτήσεις σχήματος της FEM (πολυώνυμα Lagrange). Για την μοντελοποίηση του ιδιάζοντος πεδίου τάσεων κοντά στην αιχμή της ρωγμής, εμπλουτίζουμε με τέσσερις ασυμπτωτικές συναρτήσεις, οι οποίες εκφράζονται σε ένα τοπικό πολικό σύστημα γύρω από την αιχμή και προέρχονται από αναλυτικές λύσεις της θεωρίας γραμμικής ελαστικής θραυστομηχανικής για το πεδίο μετατοπίσεων. Παράλληλα με τις εμπλουτισμένες συναρτήσεις βασης, χρειάζονται και αντίστοιχοι βαθμοί ελευθερίας σε κάθε κόμβο και για κάθε συνάρτηση εμπλουτισμού, οπότε το μέγεθος του προβλήματος αυξάνεται. Ωστόσο ο εμπλουτισμός είναι τοπικός οπότε οι επιπλέοντες βαθμοί ελευθερίας περιορίζονται γύρω από τη ρωγμή.  

Για έναν φορέα που περιέχει ασυνέχειες, η μετάβαση από την ισχυρή στην ασθενή μορφή του προβλήματος γίνεται με το θεώρημα απόκλισης, όπως στην παραδοσιακή FEM. Το προσεγγιστικό πεδίο μετατοπίσεων που αναζητείται αποτελείται από το συμβατικό πολυωνυμικό μέρος και δύο εμπλουτισμένα μέρη λόγω της Heaviside και των ασυμτπωτικών συναρτήσεων στην αιχμή. Τα μητρώα παραμορφωσης και δυσκαμψίας που προκύπτουν κατά την διακριτοποίηση της ασθενούς μορφής έχουν παρόμοια μορφή με τα αντίστοιχα μητρώα της FEM. Η διαφορά είναι ότι αποτελούνται από ξεχωριστά υπομητρώα: συμβατικά και εμπλουτισμένα. Στην XFEM δε μπορεί να χρησιμοποιηθεί η παραδοσιακή αριθμητική ολοκλήρωση Gauss, αφού οι όροι των ολοκληρωμάτων δεν είναι πια πολυωνυμικοί. Ο εμπλουτισμός με συναρτήση Heaviside εισάγει ασυνεχείς όρους στα ολοκληρώματα, ενώ ο εμπλουτισμός με τις ασυμπτωτικές συναρτήσεις στην αιχμή εισάγει ιδιάζοντες όρους. Η ολοκλήρωση Gauss είναι ακριβής μόνο για πολυωνυμικές συναρτήσεις και άρα ακατάλληλη, ακόμα και αν αυξήσουμε την τάξη της. Αντίθετα η XFEM χρησιμοποιεί δύο εναλλακτικούς κανόνες ολοκλήρωσης. Ο πρώτος βασίζεται στον χωρισμό των στοιχείων που αλληλεπιδρούν με τη ρωγμή σε υποτρίγωνα, τα οποία δεν τέμνονται από τη ρωγμή και ολοκληρώνονται με παραδοσιακού κανόνες Gauss. Έτσι η ολοκλήρωση είναι ακριβής αν το στοιχείο περιέχει μόνο πολυωνυμικούς και συνεχείς όρους. Για το δεύτερο κανόνα ολοκήρωσης, τα στοιχεία χωρίζονται σε υποτετράπλευρα που όμως μπορούν να τέμνονται από τη ρωγμή, οπότε το σφάλμα απλά περιορίζεται σε μικρότερα εμβαδά. Η μέθοδος αυτή είναι λιγότερο ακριβής, αλλά ευκολότερη στην υλοποίηση, αφού δεν απαιτεί γεωμετρικές εργασίες για την τομή της ρωγμής με το στοιχείο. Και στις δύο προσεγγίσεις η ακρίβεια της ολοκλήρωσης μπορεί να αυξηθεί χρησιμοποιώντας περισσότερα υποτρίγωνα ή υποτετράπλευρα και αυξάνοντας το σημεία ολοκλήρωσης στο εσωτιρικό αυτών.

 Η διάδοση της ρωγμής στηρίζεται στη γραμμική ελαστική θραυστομηχανική (Linear Elastic Fracture Mechanics – LEFM). Κεντρικό σημείο της LEFM είναι ο υπολογισμός των συντελεστών έντασης τάσεων. Ο συντελεστής έντασης τάσεων είναι μια ποσότητα που εκφράζει την ένταση του πεδίου τάσων γύρω από την αιχμή της ρωγμής. Ορίστηκε από τον Irwin ως το όριο του ιδιάζοντος τασικού πεδίου πλησιάζοντας την αιχμή και συνδέεται με τον ρυθμό απελευθέρωσης ενέργειας παραμόρφωσης του Griffith. Υπάρχουν τρεις δυνατές μορφές θραύσης (άνοιγμα, ολίσθηση, αποσχισμός), καθεμία πό τις οποίες έχει τον δικό της συντελεστή ένταση τάσεων. Στην πράξη εμφανίζεται μικτή μορφή θράυσης, που είναι επαλληλία των τριών βασικών μόρφων, με τους συντελεστές ένταση τάσεων να αποτελούν τους συντελεστές του γραμμικού συνδυασμού. Η παρούσα μεταπτυχιακή εργασία ασχολείται μόνο με διδιάστατα προβλήματα και δύο μορφές θραύσης.  

Οι συντελεστές έντασης τάσεων είναι απαραίτητοι για να προβλέψουμε την πορεία της ρωγμής, η οποία χαρακτηρίζεται από μια γωνία διάδοσης και ένα μήκος διάδοσης. Η γωνία διάδοσης μπορεί να υπολογιστεί με διάφορα κριτήρια, από τα οποία το πιο διαδεδομένο είναι το κριτήριο μέγιστης εφαπτομενικής τάσης. Σύμφωνα με αυτό η ρωγμή θα αναπτυχθεί σε εκείνο το επίπεδο που η εφαπτομενική τάση θα είναι κύρια, άρα μέγιστη, και η διατμητική τάση 0. Από την τελευταία σχέση προκύπτει η γωνία διάδοσης, που τελικά εξαρτάται μόνο από τον λόγο των συντελεστών έντασης τάσεων για τις δύο καταστάσεις. Το μήκος διάδοσης επιλέγεται από τον χρήστη, έτσι ώστε να είναι το ελάχιστο δυνατό, αλλά μεγαλύτερο από το μέγεθος των πεπερασμένων στοιχείων, και διατηρείται σταθερό. Αν όμοως οι συντελεστές έντασης τάσεων ξεπεράσουν μια σταθερά του υλικού που ονομάζεται αντοχή θραύσης, η ρωγμή θα διαδοθεί ακαριαία και απρόβλεπτα σε όλο το φορέα επιφέροντας κατάρρευση.  

Ο υπολογισμός των συντελεστών αντοχής τάσεων γίνεται με τη μέθοδο J-integral των Eshelby και Rice. Σύμφωνα με αυτή ο ρυθμός απελευθέρωσης ενέργειας παρομόρφωσης ισούται με ένα επικαμπύλιο ολοκλήρωμα γύρω από την αιχμή. Ο συντελεστής έντασης τάσεων για κάποια μορφή θραύσης μπορεί να υπολογιστοεί από τον ρυθμό απελευθέρωσης ενέργειας, αν υπολογίσουμε το J-integral για την επαλληλία της πραγματικης κατάστασης και μίας δυνατής κατάστασης που ταυτίζεται με τη ζητούμενη μορφή θραύσης. Το σημαντικότερο πλεονέκτημα της μεθόδου J-integral είναι ότι τα ολοκληρώματα που προκύπτουν είναι ανεξάρτητα από την καμπύλη κατά μήκος της οποίας υπολογίζονται. Για την ομαλότερη υλοποίηση της μεθόδου, τα επικαμπύλια ολοκληρώματα μετατρέπονται στη χωρική τους μορφή με τη βοήθεια του θεωρήματος απόκλισης και μιας ομαλής συνάρτησης βάρους.  

Για την υλοποίηση των παραπάνω, είναι απαραίτητος ένας τρόπος αναπαράστασης της γεωμετρίας της ρωγμής, που να εξυπηρετεί τις διάφορες γεωμετρικές ανάγκες των μεθόδων XFEM και J-integral. Μία απλή προσέγγιση είναι να περιγραφεί η ρωγμή ως μία τεθλασμένη γραμμή, στην οποία προστίθεται ένα νέο σημείο και ένα ευθύγραμμο τμήμα κάθε φορά που διαδίδεται η ρωγμή. Αυτή η μέθοδος είναι απόλυτα ακριβής αλλά οδηγεί σε χρονοβόρους αλγόριθμους, αφού για παράδειγμα πρέπει να προσπελαστούν όλα τα ευθύγραμμα τμήματα, μέχρι να βρεθεί εκείνο που τέμνει ένα πεπερασμένο στοιχείο. Επίσης χρειάζονται αλγόριθμοι υπολογιστικής γεωμετρία που δεν μπορούν να επεκταθούν εύκολα στις τρεις διαστάσεις.  

Μια εναλλακτική προσέγγιση δίνεται από τη μέθοδο συνόλων επιπέδων (Level Set Method – LSM), που αναπτύχθηκε από τους Osher & Sethian και τροποποιήθηκε για αναπαράσταση ρωγμών από τη Stolarska. Σύμφωνα με την LSM, η ρωγμή ταυτίζεται με τη μηδενική ισοϋψή μιας συνάρτησης φ που επιστρέφει την προσημασμένη απόσταση ενός σημείου από τη ρωγμή. Η αιχμή της ρωγμής προσδιορίζεται χρησιμοποιώντας μια δεύτερη συνάρτηση ψ, η μηδενική ισοϋψής της οποίας είναι κάθετη στην φ=0 στο σημείο της αιχμής. Αυτές οι συναρτήσεις αποθηκεύονται στους κόμβους του πλέγματος και παρεμβάλλονται στο εσωτερικό των πεπερασμένων στοιχείων με τις συναρτήσεις σχήματος. Έτσι προκύπτουν κομψοί και πολύ αποδοτικοί αλγόριθμου, όταν συνδυάζεται με την XFEM. Ταυτόχρονα η LSM μπορεί εύκολα να επεκταθεί σε τρισδιάστατα προβλήματα. Ωστόσο δεν είναι απόλυτα ακριβής, αφού δεν μπορεί να αναπαραστήσει την αλλαγή κλίσης της ρωγμής στο εσωτερικό των στοιχείων, λόγω του ότι στηρίζεται στις γραμμικές συναρτήσεις σχήματος.  

Όλα τα παραπάνω έχουν υλοποιηθεί για διδιάστατα προβλήματα στην αντικειμενοστραφή γλώσσα προγραμματισμού C#. Ο κώδικας αποτελεί μέρος του λογισμικού MSolve που αναπτύσσεται από το Εργαστήριο Στατικής και Αντισεισμικών Ερευνών, στη σχολή Πολιτικών Μηχανικών, ΕΜΠ. Τέλος έχουν εφαρμοστεί σε μια σειρά από αριθμητικά παραδείγματα, από τα οποία προκύπτουν δύο βασικά συμπεράσματα. Πρώτον, ότι η μέθοδος XFEM είναι δυνατή να δώσης λύσεις υψηλής ακρίβειας που ταυτίζονται με αναλυτικές επιλύσεις, όπου αυτές είναι διαθέσιμες, και δεύτερον, ότι η επιλογή του μήκους διάδοσης της ρωγμής πρέπει να συνεκτιμάται με την επιλογή του μεγέθους των πεπερασμένων στοιχείων και της ακτίνας του J-integral. 

 

Δείτε τη ΜΕ στη βιβλιοθήκη του ΕΜΠ