Перейти до змісту

Метод Ньютона—Рафсона

Лінеаризація та ітераційна рекурентна формула

Нелінійне рівняння для вузлового переміщення \(\boldsymbol{u}_{n+1}\) у момент часу \(t_{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}}\) — представлення тензора визначального співвідношення \(\boldsymbol{\mathsf{C}}\) у позначеннях Voigt, тобто матриця матеріальної жорсткості \(6\times 6\) (Тензорні позначення та математичні основи). \(\boldsymbol{S}_9, \boldsymbol{F}_9\) — наведені нижче матриці перестановки, що використовуються для подання члена геометричної жорсткості як матричного добутку. Спочатку для тензора другого рангу \(\boldsymbol{A}\) розміру \(3\times 3\) визначимо позначення \([\,\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\) на блоки \(\boldsymbol{K}^e_{\alpha\beta}\) розміру \(d\times d\) для кожної пари вузлів і використанням множини складання тензорів другого рангу \(\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