Иди на текст

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

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

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

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

AI-assisted translation May contain errors Official docs Status