לדלג לתוכן

שיטת ניוטון–רפסון

לינאריזציה ונוסחת הנסיגה האיטרטיבית

את המשוואה הלא־לינארית עבור ההזזה הצומתית בזמן \(t_{n+1}\),‏ \(\boldsymbol{u}_{n+1}\), שהתקבלה בעבודה הווירטואלית של הכוחות החיצוניים והרכבת המשוואה הגלובלית, פותרים בשיטת ניוטון–רפסון. מניחים שעד לזמן \(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}\), ומתקדמים לצעד הזמן הבא.

נושאים קשורים