Μετάβαση στο περιεχόμενο

Μέθοδος Newton-Raphson

Γραμμικοποίηση και επαναληπτική αναδρομή

Η Εικονική εργασία εξωτερικών δυνάμεων και συναρμολόγηση της καθολικής εξίσωσης δίνει μια μη γραμμική εξίσωση για την κομβική μετατόπιση στη χρονική στιγμή \(t_{n+1}\), \(\boldsymbol{u}_{n+1}\), η οποία επιλύεται με τη μέθοδο Newton-Raphson. Η κομβική μετατόπιση έως τη χρονική στιγμή \(t_n\), \(\boldsymbol{u}_n\), θεωρείται γνωστή και η αύξηση μετατόπισης \(\Delta\boldsymbol{u}\) λαμβάνεται ως η άγνωστη μεταβλητή που πρέπει να προσδιοριστεί

\[ \boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u} \]

Στη συνέχεια παραλείπεται η εξάρτηση του διανύσματος εξωτερικής δύναμης από την κομβική μετατόπιση και, με \(\boldsymbol{F}(\boldsymbol{u}_{n+1}) = \boldsymbol{F}_{n+1}\),

\[ \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u}) - \boldsymbol{F}_{n+1} = \boldsymbol{0} \]

επιλύεται.

Στην τρέχουσα λύση \(\Delta\boldsymbol{u}\), ορίζεται η εφαπτομενική δυσκαμψία

\[ \boldsymbol{K} = \left. \frac{\partial \boldsymbol{Q}}{\partial \boldsymbol{u}} \right|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}} \]

Με αυτήν, η γραμμικοποίηση της μη γραμμικής εξίσωσης δίνει

\[ \boldsymbol{K}\, d\boldsymbol{u} + \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u}) - \boldsymbol{F}_{n+1} = \boldsymbol{0} \]

Έστω ότι η διόρθωση στην \(i\)-οστή επανάληψη είναι \(d\boldsymbol{u}_i\) και ότι το διάνυσμα υπολοίπου στην αρχή της επανάληψης είναι

\[ \boldsymbol{R}_{i-1} = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u}) \]

Τότε η επαναληπτική αναδρομή είναι

\[ \boldsymbol{K}_i\, d\boldsymbol{u}_i = \boldsymbol{R}_{i-1}, \qquad \Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i \]

Επομένως, το υπόλοιπο \(\boldsymbol{R}_i\) είναι ποσότητα που αντιστοιχεί στην ανισορροπία δυνάμεων ως προς την ισορροπία.

Κατασκευή του πίνακα εφαπτομενικής δυσκαμψίας

Η εφαπτομενική δυσκαμψία \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\) κατασκευάζεται με μερική παραγώγιση του διανύσματος εσωτερικών δυνάμεων του στοιχείου που προκύπτει στην Διακριτοποίηση της εικονικής εργασίας των εσωτερικών δυνάμεων ως προς την κομβική μετατόπιση, ολοκλήρωση των ολοκληρωτέων σε επίπεδο στοιχείου πάνω σε κάθε περιοχή στοιχείου και συναρμολόγησή τους. Συμβολίζοντας το ολοκληρωτέο σε επίπεδο στοιχείου με \(\boldsymbol{K}^e_X\) (συμβολισμός διαμόρφωσης αναφοράς, διατύπωση TL) ή \(\boldsymbol{K}^e_x\) (συμβολισμός τρέχουσας διαμόρφωσης, διατύπωση UL), η εφαπτομενική δυσκαμψία του στοιχείου είναι

\[ \boldsymbol{K}^e = \int_{\Omega^e_0} \boldsymbol{K}^e_X\, dV \quad (\text{TL}), \qquad \boldsymbol{K}^e = \int_{\Omega^e} \boldsymbol{K}^e_x\, dv \quad (\text{UL}) \]

Παρακάτω δίνονται οι τελικές μορφές των ολοκληρωτέων TL/UL. Και στις δύο περιπτώσεις αναλύονται στο άθροισμα ενός όρου δυσκαμψίας υλικού (όρος αρχικής μετατόπισης) και ενός όρου γεωμετρικής δυσκαμψίας (όρος αρχικής τάσης).

Διατύπωση Total Lagrange

Στη διατύπωση Total Lagrange θεωρείται γραμμική σχέση μεταξύ του ρυθμού της δεύτερης τάσης Piola-Kirchhoff \(\dot{\boldsymbol{S}}\) και του ρυθμού παραμόρφωσης Green-Lagrange \(\dot{\boldsymbol{E}}\), δηλαδή \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\). Αυτό αντιστοιχεί σε καταστατικούς νόμους για γραμμικά ελαστικά υλικά (υλικά St. Venant-Kirchhoff) και υπερελαστικά υλικά, και το FrontISTR χρησιμοποιεί τη διατύπωση Total Lagrange για αυτά τα υλικά. Το ολοκληρωτέο της εφαπτομενικής δυσκαμψίας του στοιχείου γράφεται τότε σε τανυστική μορφή ως

\[ \delta\boldsymbol{u}^{eT}\, \boldsymbol{K}^e_X\, \dot{\boldsymbol{u}}^e = \dot{\boldsymbol{S}}:\delta\boldsymbol{E} + \boldsymbol{S}:(\delta\boldsymbol{F}^T \dot{\boldsymbol{F}}) \]

Ο πρώτος όρος στο δεξιό μέλος είναι ο όρος δυσκαμψίας υλικού (όρος αρχικής μετατόπισης) και ο δεύτερος είναι ο όρος γεωμετρικής δυσκαμψίας (όρος αρχικής τάσης).

Στην υλοποίηση του FrontISTR, αυτό το ολοκληρωτέο υπολογίζεται σε μορφή πίνακα με χρήση του συμβολισμού Voigt:

\[ \boldsymbol{K}^e_X = (\boldsymbol{B}_L + \boldsymbol{B}_{NL})^T\, \tilde{\boldsymbol{C}}\, (\boldsymbol{B}_L + \boldsymbol{B}_{NL}) + \boldsymbol{F}_9^T\, \boldsymbol{S}_9\, \boldsymbol{F}_9 \]

Οι πίνακες είναι οι εξής. Τα \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}\) είναι οι πίνακες B που εισάγονται στη Διακριτοποίηση της εικονικής εργασίας των εσωτερικών δυνάμεων, και το \(\tilde{\boldsymbol{C}}\) είναι η αναπαράσταση Voigt του καταστατικού τανυστή \(\boldsymbol{\mathsf{C}}\), δηλαδή ένας πίνακας δυσκαμψίας υλικού \(6\times 6\) (Τανυστικός συμβολισμός και μαθηματικά θεμέλια). Τα \(\boldsymbol{S}_9, \boldsymbol{F}_9\) είναι οι ακόλουθοι πίνακες αναδιάταξης που χρησιμοποιούνται για την έκφραση του όρου γεωμετρικής δυσκαμψίας ως γινομένου πινάκων. Αρχικά, για έναν τανυστή δεύτερης τάξης \(3\times 3\), \(\boldsymbol{A}\), ορίζεται ο συμβολισμός \([\,\cdot\,]\) που τον αναδιατάσσει σε διάνυσμα 9 συνιστωσών ως

\[ [\boldsymbol{A}] = (A_{11}, A_{21}, A_{31}, A_{12}, A_{22}, A_{32}, A_{13}, A_{23}, A_{33})^T \]

Με αυτόν τον ορισμό, το \(\boldsymbol{F}_9\) εκφράζει τη μεταβολή της βαθμίδας παραμόρφωσης με τη μορφή \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\) και είναι πίνακας \(9\times d n_e\). Για τον κόμβο στοιχείου \(\alpha = 1, \ldots, n_e\), το αντίστοιχο μπλοκ \(9\times d\) είναι

\[ [\boldsymbol{F}_9]_\alpha = \begin{bmatrix} (\partial N_\alpha^e/\partial X_1)\, \boldsymbol{I} \\ (\partial N_\alpha^e/\partial X_2)\, \boldsymbol{I} \\ (\partial N_\alpha^e/\partial X_3)\, \boldsymbol{I} \end{bmatrix} \qquad (\boldsymbol{I} \text{ είναι ο } 3\times 3 \text{ μοναδιαίος πίνακας}) \]

και δίνεται από \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\), με τα μπλοκ τοποθετημένα οριζόντια με τη σειρά των κόμβων του στοιχείου. Το \(\boldsymbol{S}_9\) επιλέγεται έτσι ώστε, σε συνδυασμό με αυτόν τον πίνακα, ο όρος γεωμετρικής δυσκαμψίας να εκφράζεται ως \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\)· είναι ο ακόλουθος πίνακας \(9\times 9\)

\[ \boldsymbol{S}_9 = \begin{bmatrix} S_{11} \boldsymbol{I} & S_{12} \boldsymbol{I} & S_{13} \boldsymbol{I} \\ S_{21} \boldsymbol{I} & S_{22} \boldsymbol{I} & S_{23} \boldsymbol{I} \\ S_{31} \boldsymbol{I} & S_{32} \boldsymbol{I} & S_{33} \boldsymbol{I} \end{bmatrix} \]

Αυτός είναι ο προκύπτων πίνακας.

Διατύπωση Updated Lagrange

Στη διατύπωση Updated Lagrange θεωρείται γραμμική σχέση μεταξύ του ρυθμού Jaumann του σχετικού τανυστή τάσης Kirchhoff \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) και του τανυστή ρυθμού παραμόρφωσης \(\boldsymbol{D}\), δηλαδή \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\). Αυτή είναι η μορφή ενός υποελαστικού καταστατικού νόμου κοινή για γραμμικά ελαστικά, ελαστοπλαστικά υλικά και υλικά ερπυσμού, και το FrontISTR χρησιμοποιεί τη διατύπωση Updated Lagrange για αυτά τα υλικά. Το ολοκληρωτέο της εφαπτομενικής δυσκαμψίας του στοιχείου, εκφρασμένο στην τρέχουσα διαμόρφωση, γράφεται τότε σε τανυστική μορφή ως

\[ \delta\boldsymbol{u}^{eT}\, \boldsymbol{K}^e_x\, \dot{\boldsymbol{u}}^e = \boldsymbol{\sigma}^{\nabla T}:\delta\boldsymbol{A}_{(L)} + \boldsymbol{\sigma}:(\delta\boldsymbol{F}_t^T\, \boldsymbol{L}) \]

όπου \(\boldsymbol{\sigma}^{\nabla T}\) είναι ο ρυθμός Truesdell, \(\boldsymbol{A}_{(L)}\) είναι το γραμμικό μέρος της παραμόρφωσης Almansi, \(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\) είναι η βαθμίδα μετατόπισης ως προς την τρέχουσα διαμόρφωση και \(\boldsymbol{L}\) είναι ο τανυστής βαθμίδας ταχύτητας. Ο πρώτος όρος στο δεξιό μέλος είναι ο όρος δυσκαμψίας υλικού και ο δεύτερος είναι ο όρος γεωμετρικής δυσκαμψίας.

Στην υλοποίηση του FrontISTR, αυτό το ολοκληρωτέο υπολογίζεται σε μορφή πίνακα με χρήση του συμβολισμού Voigt:

\[ \boldsymbol{K}^e_x = \boldsymbol{b}^T\, (\tilde{\boldsymbol{C}} - \boldsymbol{G})\, \boldsymbol{b} + \boldsymbol{f}_9^T\, \boldsymbol{\sigma}_9\, \boldsymbol{f}_9 \]

Εδώ, το \(\boldsymbol{b}\) είναι ο πίνακας B που κατασκευάζεται στην τρέχουσα διαμόρφωση (Διακριτοποίηση της εικονικής εργασίας των εσωτερικών δυνάμεων). Τα \(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\) προκύπτουν από τα \(\boldsymbol{S}_9, \boldsymbol{F}_9\) που ορίστηκαν για τη διατύπωση TL, αντικαθιστώντας τη δεύτερη τάση PK \(\boldsymbol{S}\) με την τάση Cauchy \(\boldsymbol{\sigma}\) και τη βαθμίδα διαμόρφωσης αναφοράς \(\partial N_\alpha^e/\partial X_i\) με τη βαθμίδα τρέχουσας διαμόρφωσης \(\partial N_\alpha^e/\partial x_i\).

Το \(\boldsymbol{G}\) είναι πίνακας διόρθωσης εξαρτώμενος από την τάση Cauchy, ο οποίος απαιτείται ώστε ο υποελαστικός καταστατικός νόμος \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) να είναι συνεπής με το πλαίσιο εφαπτομενικής δυσκαμψίας ως καταστατικός νόμος βασισμένος στον ρυθμό Truesdell. Προκύπτει με διάταξη των συνιστωσών του τανυστή τέταρτης τάξης \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) σε μορφή Voigt \(6\times 6\) ως

\[ \boldsymbol{G} = \begin{bmatrix} 2\sigma_{11} & 0 & 0 & \sigma_{12} & 0 & \sigma_{31} \\ 0 & 2\sigma_{22} & 0 & \sigma_{12} & \sigma_{23} & 0 \\ 0 & 0 & 2\sigma_{33} & 0 & \sigma_{23} & \sigma_{31} \\ \sigma_{12} & \sigma_{12} & 0 & \tfrac{\sigma_{11}+\sigma_{22}}{2} & \tfrac{\sigma_{12}}{2} & \tfrac{\sigma_{23}}{2} \\ 0 & \sigma_{23} & \sigma_{23} & \tfrac{\sigma_{12}}{2} & \tfrac{\sigma_{22}+\sigma_{33}}{2} & \tfrac{\sigma_{12}}{2} \\ \sigma_{31} & 0 & \sigma_{31} & \tfrac{\sigma_{23}}{2} & \tfrac{\sigma_{12}}{2} & \tfrac{\sigma_{33}+\sigma_{11}}{2} \end{bmatrix} \]

Αυτός είναι ο προκύπτων πίνακας.

Συναρμολόγηση του καθολικού πίνακα δυσκαμψίας

Η καθολική εφαπτομενική δυσκαμψία \(\boldsymbol{K}\) προκύπτει διαιρώντας κάθε δυσκαμψία στοιχείου \(\boldsymbol{K}^e\) σε μπλοκ \(d\times d\), \(\boldsymbol{K}^e_{\alpha\beta}\), για κάθε ζεύγος κόμβων και χρησιμοποιώντας το σύνολο συναρμολόγησης τανυστών δεύτερης τάξης \(\mathcal{E}^2(i_g, i_h)\) που εισάγεται στη Συναρμολόγηση φυσικών μεγεθών κόμβων στοιχείου:

\[ \boldsymbol{K}_{i_gi_h} = \sum_{(e,\alpha,\beta) \in \mathcal{E}^2(i_g, i_h)} \boldsymbol{K}^e_{\alpha\beta} \]

Οι προκύπτουσες τιμές διατάσσονται σε πίνακα με γραμμή \(i_g\) και στήλη \(i_h\). Στην υλοποίηση, το σύνολο \(\mathcal{E}^2\) δεν κατασκευάζεται ρητά· αντί αυτού, τα αντίστοιχα μπλοκ προστίθενται απευθείας μέσα στον βρόχο στοιχείων. Ο πίνακας είναι τετραγωνικός με διάσταση ίση με βαθμούς ελευθερίας ανά κόμβο \(\times\) συνολικό αριθμό κόμβων \(n_g\), αλλά επειδή οι συνιστώσες εκτός εκείνων μεταξύ κόμβων που συνδέονται μέσω στοιχείων είναι \(0\), αποθηκεύεται σε αραιή μορφή.

Οι πίνακες δυσκαμψίας στοιχείων για τις διατυπώσεις TL και UL έχουν την ίδια μορφή, με εξαίρεση την αλλαγή της διαμόρφωσης αναφοράς (κομβικές συντεταγμένες και πηγή που χρησιμοποιείται για την κατασκευή του πίνακα B) και την παρουσία ή απουσία του πίνακα \(\boldsymbol{G}\). Επομένως, το FrontISTR υλοποιεί και τις δύο διατυπώσεις σε κοινή υπορουτίνα.

Αλγόριθμος επανάληψης

Συνοψίζοντας τα παραπάνω, στην αρχή της επανάληψης θέστε \(\Delta\boldsymbol{u} = \boldsymbol{0}\) και υπολογίστε το αρχικό υπόλοιπο \(\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)\). Στη συνέχεια, στην \(i\)-οστή επανάληψη, εκτελέστε την ακόλουθη διαδικασία.

  1. Στην τρέχουσα μετατόπιση \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\), υπολογίστε την εφαπτομενική δυσκαμψία \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\) με τη διαδικασία της Κατασκευής του πίνακα εφαπτομενικής δυσκαμψίας.
  2. Για την επιβολή γεωμετρικών οριακών συνθηκών, τροποποιήστε τον πίνακα εφαπτομενικής δυσκαμψίας και το διάνυσμα υπολοίπου για τους βαθμούς ελευθερίας που υπόκεινται σε περιορισμούς μετατόπισης, λαμβάνοντας \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\) (Χειρισμός γεωμετρικών οριακών συνθηκών).
  3. Επιλύστε τη γραμμική εξίσωση \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\) για να λάβετε τη διόρθωση \(d\boldsymbol{u}_i\). Η διαδικασία αυτή συχνά αντιπροσωπεύει το μεγαλύτερο μέρος του υπολογιστικού κόστους της επαναληπτικής επίλυσης.
  4. Ενημερώστε την αύξηση μετατόπισης ως \(\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i\) και, αντίστοιχα, υπολογίστε το διάνυσμα εσωτερικών δυνάμεων \(\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) και το υπόλοιπο \(\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\).
  5. Ελέγξτε τη σύγκλιση και τερματίστε την επανάληψη εάν επιτευχθεί σύγκλιση. Στο υπόλοιπο \(\boldsymbol{R}_i\) εμφανίζονται συνιστώσες που αντιστοιχούν σε αντιδράσεις περιορισμών στους βαθμούς ελευθερίας που υπόκεινται σε γεωμετρικές οριακές συνθήκες, επομένως ο δείκτης σύγκλισης κατασκευάζεται από το \(\tilde{\boldsymbol{R}}_i\) αφού εξαιρεθούν αυτές οι συνιστώσες. Συγκεκριμένοι δείκτες και κατώφλια σύγκλισης περιγράφονται στα Κριτήρια σύγκλισης. Εάν δεν επιτευχθεί σύγκλιση και επιτευχθεί το όριο επαναλήψεων, η επανάληψη θεωρείται αποτυχημένη.

Όταν η επανάληψη συγκλίνει, προσθέστε το συγκλίνον \(\Delta\boldsymbol{u}\) στο \(\boldsymbol{u}_n\) για να λάβετε τη συσσωρευμένη μετατόπιση στη χρονική στιγμή \(t_{n+1}\), \(\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}\), και προχωρήστε στο επόμενο χρονικό βήμα.

Σχετικά θέματα

AI-assisted translation May contain errors Official docs Status