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

Њутн-Рафсонов метод

Линеаризација и итеративна рекурзија

Виртуелна работа на надворешните сили и составување на глобалната равенка дава нелинеарна равенка за јазолното поместување во времето \(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}}\). Ова одговара на конститутивни закони за линеарно еластични материјали (материјали 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, овој подинтегрален израз се пресметува во матрична форма со Војтова нотација:

\[ \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 \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}\), и продолжете на следниот временски чекор.

Поврзани теми