پرش به محتویات

روش Newton-Raphson

خطی‌سازی و رابطه بازگشتی تکرار

معادله غیرخطی مربوط به جابه‌جایی گره‌ای \(\boldsymbol{u}_{n+1}\) در زمان \(t_{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 ارائه می‌شود. در هر دو حالت، عبارت به مجموع جمله سختی ماده (جمله جابه‌جایی اولیه) و جمله سختی هندسی (جمله تنش اولیه) تجزیه می‌شود.

روش لاگرانژ کل

در روش لاگرانژ کل، رابطه خطی \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\) میان نرخ تنش Piola-Kirchhoff دوم \(\dot{\boldsymbol{S}}\) و نرخ کرنش Green-Lagrange \(\dot{\boldsymbol{E}}\) فرض می‌شود. این رابطه متناظر با قوانین رفتاری مواد کشسان خطی (St. Venant-Kirchhoff) و مواد هایپرالاستیک است و FrontISTR روش لاگرانژ کل را برای این مواد به‌کار می‌برد. در این حالت، عبارت زیرانتگرالی سختی مماسی المان به‌صورت تانسوری چنین نوشته می‌شود:

\[ \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\) ماتریس \(9\times 9\) زیر است که به‌گونه‌ای انتخاب شده است تا همراه با این ماتریس، جمله سختی هندسی به‌صورت \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\) بیان شود:

\[ \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} \]

است.

روش لاگرانژ به‌روزشده

در روش لاگرانژ به‌روزشده، رابطه خطی \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) میان نرخ Jaumann تانسور تنش Kirchhoff نسبی \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) و تانسور نرخ تغییرشکل \(\boldsymbol{D}\) فرض می‌شود. این شکل قانون رفتاری hypoelastic در مواد کشسان خطی، الاستوپلاستیک و خزشی مشترک است و FrontISTR روش لاگرانژ به‌روزشده را برای این مواد به‌کار می‌برد. در این حالت، عبارت زیرانتگرالی سختی مماسی المان در پیکربندی جاری به‌صورت تانسوری چنین نوشته می‌شود:

\[ \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 است که برای سازگار کردن قانون رفتاری hypoelastic \(\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}\)، به‌دست آید و سپس تحلیل به گام زمانی بعدی می‌رود.

موارد مرتبط