Перейти к содержанию

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

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

Нелинейное уравнение в момент времени \(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. В обоих случаях они раскладываются в сумму материального члена жёсткости (члена начального перемещения) и геометрического члена жёсткости (члена начального напряжения).

Полная лагранжева формулировка

В полной лагранжевой формулировке предполагается линейная связь между скоростью второго напряжения Пиолы–Кирхгофа \(\dot{\boldsymbol{S}}\) и скоростью деформации Грина–Лагранжа \(\dot{\boldsymbol{E}}\): \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\). Это соответствует определяющим законам для линейно-упругих материалов (материалов Сен-Венана–Кирхгофа) и гиперупругих материалов; в 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 это подынтегральное выражение вычисляется в матричной форме с использованием записи Фойгта:

\[ \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}}\) в записи Фойгта, то есть матрица материальной жёсткости \(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} \]

Это полученная матрица.

Обновлённая лагранжева формулировка

В обновлённой лагранжевой формулировке предполагается линейная связь между производной Яуманна относительного тензора напряжений Кирхгофа \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) и тензором скоростей деформации \(\boldsymbol{D}\): \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\). Это форма гипоупругого определяющего закона, общая для линейно-упругих, упругопластических и ползучих материалов; в 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}\) — производная Трусделла, \(\boldsymbol{A}_{(L)}\) — линейная часть деформации Альманси, \(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\) — градиент перемещения относительно текущей конфигурации, а \(\boldsymbol{L}\) — тензор градиента скорости. Первый член правой части — материальный член жёсткости, второй — геометрический член жёсткости.

В реализации FrontISTR это подынтегральное выражение вычисляется в матричной форме с использованием записи Фойгта:

\[ \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}\) на напряжение Коши \(\boldsymbol{\sigma}\) и градиента по исходной конфигурации \(\partial N_\alpha^e/\partial X_i\) на градиент по текущей конфигурации \(\partial N_\alpha^e/\partial x_i\).

\(\boldsymbol{G}\) — зависящая от напряжения Коши поправочная матрица, необходимая для согласования гипоупругого определяющего закона \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) с формулировкой касательной жёсткости как определяющего закона на основе производной Трусделла. Она получается представлением компонент тензора четвёртого ранга \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) в записи Фойгта \(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}\), после чего выполняется переход к следующему временному шагу.

Связанные темы