Към съдържанието

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

Линеаризация и итерационна рекурентност

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

\[ \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\) са следните матрици за пренареждане, използвани за представяне на члена на геометричната коравина като матрично произведение. Първо, за тензор от втори ранг \(\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) \]

и \(\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 се приема линейна зависимост между скоростта на Яуман на тензора на относителното напрежение на Кирхоф \(\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}\) е скоростта на Трюсдел, \(\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}\), след което се преминава към следващата времева стъпка.

Свързани теми

AI-assisted translation May contain errors Official docs Status