Przejdź do treści

Metoda Newtona-Raphsona

Linearyzacja i rekurencja iteracyjna

Praca wirtualna sił zewnętrznych i składanie równania globalnego prowadzi do nieliniowego równania względem przemieszczenia węzłowego w chwili \(t_{n+1}\), \(\boldsymbol{u}_{n+1}\), które jest rozwiązywane metodą Newtona-Raphsona. Przyjmuje się, że przemieszczenie węzłowe do chwili \(t_n\), \(\boldsymbol{u}_n\), jest znane, a przyrost przemieszczenia \(\Delta\boldsymbol{u}\) stanowi niewiadomą do wyznaczenia

\[ \boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u} \]

Dalej pomija się zależność wektora sił zewnętrznych od przemieszczenia węzłowego i przy \(\boldsymbol{F}(\boldsymbol{u}_{n+1}) = \boldsymbol{F}_{n+1}\),

\[ \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u}) - \boldsymbol{F}_{n+1} = \boldsymbol{0} \]

rozwiązuje się.

Dla bieżącego rozwiązania \(\Delta\boldsymbol{u}\) zdefiniuj sztywność styczną

\[ \boldsymbol{K} = \left. \frac{\partial \boldsymbol{Q}}{\partial \boldsymbol{u}} \right|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}} \]

Korzystając z tego, linearyzacja równania nieliniowego daje

\[ \boldsymbol{K}\, d\boldsymbol{u} + \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u}) - \boldsymbol{F}_{n+1} = \boldsymbol{0} \]

Niech poprawka w \(i\)-tej iteracji będzie równa \(d\boldsymbol{u}_i\), a wektor residuum na początku iteracji będzie równy

\[ \boldsymbol{R}_{i-1} = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u}) \]

Wtedy rekurencja iteracyjna ma postać

\[ \boldsymbol{K}_i\, d\boldsymbol{u}_i = \boldsymbol{R}_{i-1}, \qquad \Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i \]

Zatem residuum \(\boldsymbol{R}_i\) jest wielkością odpowiadającą niezrównoważeniu sił względem stanu równowagi.

Konstrukcja macierzy sztywności stycznej

Sztywność styczna \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\) jest konstruowana przez częściowe zróżniczkowanie wektora sił wewnętrznych elementu uzyskanego w Dyskretyzacji pracy wirtualnej sił wewnętrznych względem przemieszczenia węzłowego, całkowanie otrzymanych wyrażeń na poziomie elementu po domenie każdego elementu oraz ich składanie. Oznaczając wyrażenie na poziomie elementu jako \(\boldsymbol{K}^e_X\) (notacja konfiguracji odniesienia, sformułowanie TL) lub \(\boldsymbol{K}^e_x\) (notacja konfiguracji bieżącej, sformułowanie UL), sztywność styczna elementu wynosi

\[ \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}) \]

Poniżej podano ostateczne postacie wyrażeń TL/UL. W obu przypadkach są one rozkładane na sumę składnika sztywności materiałowej (składnika przemieszczenia początkowego) oraz składnika sztywności geometrycznej (składnika naprężenia początkowego).

Sformułowanie całkowite Lagrange’a

W sformułowaniu całkowitym Lagrange’a zakłada się liniową zależność między szybkością drugiego naprężenia Pioli-Kirchhoffa \(\dot{\boldsymbol{S}}\) a szybkością odkształcenia Greena-Lagrange’a \(\dot{\boldsymbol{E}}\), mianowicie \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\). Odpowiada to prawom konstytutywnym dla materiałów liniowo sprężystych (materiałów St. Venanta-Kirchhoffa) i materiałów hipersprężystych, a FrontISTR używa dla tych materiałów sformułowania całkowitego Lagrange’a. Wyrażenie sztywności stycznej elementu zapisuje się wtedy w postaci tensorowej jako

\[ \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}}) \]

Pierwszy składnik po prawej stronie jest składnikiem sztywności materiałowej (składnikiem przemieszczenia początkowego), a drugi — składnikiem sztywności geometrycznej (składnikiem naprężenia początkowego).

W implementacji FrontISTR wyrażenie to jest obliczane w postaci macierzowej przy użyciu notacji Voigta:

\[ \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 \]

Poszczególne macierze są następujące. \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}\) są macierzami B wprowadzonymi w Dyskretyzacji pracy wirtualnej sił wewnętrznych, a \(\tilde{\boldsymbol{C}}\) jest reprezentacją Voigta tensora konstytutywnego \(\boldsymbol{\mathsf{C}}\), czyli macierzą sztywności materiałowej \(6\times 6\) (Notacja tensorowa i podstawy matematyczne). \(\boldsymbol{S}_9, \boldsymbol{F}_9\) są następującymi macierzami przestawienia używanymi do wyrażenia składnika sztywności geometrycznej jako iloczynu macierzy. Najpierw, dla tensora drugiego rzędu \(3\times 3\) \(\boldsymbol{A}\), definiuje się zapis \([\,\cdot\,]\), który przestawia go do wektora 9-składowego jako

\[ [\boldsymbol{A}] = (A_{11}, A_{21}, A_{31}, A_{12}, A_{22}, A_{32}, A_{13}, A_{23}, A_{33})^T \]

Przy tej definicji \(\boldsymbol{F}_9\) wyraża wariację gradientu deformacji w postaci \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\) i jest macierzą \(9\times d n_e\). Dla węzła elementu \(\alpha = 1, \ldots, n_e\) odpowiadający blok \(9\times d\) ma postać

\[ [\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{ jest macierzą jednostkową } 3\times 3) \]

i jest dane przez \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\), przy czym bloki są ułożone poziomo w kolejności węzłów elementu. \(\boldsymbol{S}_9\) jest dobierana tak, aby w połączeniu z tą macierzą składnik sztywności geometrycznej był wyrażony jako \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\); jest to następująca macierz \(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} \]

Oto otrzymana macierz.

Sformułowanie zaktualizowane Lagrange’a

W sformułowaniu zaktualizowanym Lagrange’a zakłada się liniową zależność między pochodną Jaumanna względnego tensora naprężenia Kirchhoffa \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) a tensorem szybkości odkształcenia \(\boldsymbol{D}\), mianowicie \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\). Jest to postać hiposprężystego prawa konstytutywnego wspólna dla materiałów liniowo sprężystych, sprężysto-plastycznych i pełzających, a FrontISTR używa dla tych materiałów sformułowania zaktualizowanego Lagrange’a. Wyrażenie sztywności stycznej elementu zapisane w konfiguracji bieżącej ma wtedy w postaci tensorowej postać

\[ \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}) \]

gdzie \(\boldsymbol{\sigma}^{\nabla T}\) jest pochodną Truesdella, \(\boldsymbol{A}_{(L)}\) jest liniową częścią odkształcenia Almansiego, \(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\) jest gradientem przemieszczenia względem konfiguracji bieżącej, a \(\boldsymbol{L}\) jest tensorem gradientu prędkości. Pierwszy składnik po prawej stronie jest składnikiem sztywności materiałowej, a drugi — składnikiem sztywności geometrycznej.

W implementacji FrontISTR wyrażenie to jest obliczane w postaci macierzowej przy użyciu notacji Voigta:

\[ \boldsymbol{K}^e_x = \boldsymbol{b}^T\, (\tilde{\boldsymbol{C}} - \boldsymbol{G})\, \boldsymbol{b} + \boldsymbol{f}_9^T\, \boldsymbol{\sigma}_9\, \boldsymbol{f}_9 \]

Tutaj \(\boldsymbol{b}\) jest macierzą B skonstruowaną w konfiguracji bieżącej (Dyskretyzacja pracy wirtualnej sił wewnętrznych). \(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\) otrzymuje się z \(\boldsymbol{S}_9, \boldsymbol{F}_9\) zdefiniowanych dla sformułowania TL przez zastąpienie drugiego naprężenia PK \(\boldsymbol{S}\) naprężeniem Cauchy’ego \(\boldsymbol{\sigma}\) oraz gradientu konfiguracji odniesienia \(\partial N_\alpha^e/\partial X_i\) gradientem konfiguracji bieżącej \(\partial N_\alpha^e/\partial x_i\).

\(\boldsymbol{G}\) jest macierzą korekcyjną zależną od naprężenia Cauchy’ego, wymaganą do zapewnienia zgodności hiposprężystego prawa konstytutywnego \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) z ramami sztywności stycznej jako prawa konstytutywnego opartego na pochodnej Truesdella. Otrzymuje się ją przez uporządkowanie składowych tensora czwartego rzędu \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) w postaci Voigta \(6\times 6\) jako

\[ \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} \]

Oto otrzymana macierz.

Składanie globalnej macierzy sztywności

Globalną sztywność styczną \(\boldsymbol{K}\) otrzymuje się przez podział każdej sztywności elementu \(\boldsymbol{K}^e\) na bloki \(d\times d\) \(\boldsymbol{K}^e_{\alpha\beta}\) dla każdej pary węzłów i użycie zbioru składania tensorów drugiego rzędu \(\mathcal{E}^2(i_g, i_h)\) wprowadzonego w Składaniu fizycznych wielkości elementowo-węzłowych:

\[ \boldsymbol{K}_{i_gi_h} = \sum_{(e,\alpha,\beta) \in \mathcal{E}^2(i_g, i_h)} \boldsymbol{K}^e_{\alpha\beta} \]

Otrzymane wartości są rozmieszczane w macierzy z wierszem \(i_g\) i kolumną \(i_h\). W implementacji zbiór \(\mathcal{E}^2\) nie jest jawnie konstruowany; zamiast tego odpowiadające bloki są dodawane bezpośrednio w pętli po elementach. Macierz jest kwadratowa, o wymiarze równym liczbie stopni swobody na węzeł \(\times\) całkowitej liczbie węzłów \(n_g\), lecz ponieważ składowe inne niż między węzłami połączonymi przez elementy są równe \(0\), jest przechowywana w postaci macierzy rzadkiej.

Macierze sztywności elementów dla sformułowań TL i UL mają tę samą postać z wyjątkiem przełączenia konfiguracji odniesienia (współrzędnych węzłowych i źródła używanego do konstrukcji macierzy B) oraz obecności lub braku macierzy \(\boldsymbol{G}\). FrontISTR implementuje zatem oba sformułowania we wspólnym podprogramie.

Algorytm iteracji

Podsumowując powyższe, na początku iteracji ustaw \(\Delta\boldsymbol{u} = \boldsymbol{0}\) i oblicz początkowe residuum \(\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)\). Następnie w \(i\)-tej iteracji wykonaj następującą procedurę.

  1. Dla bieżącego przemieszczenia \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\) oblicz sztywność styczną \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\) zgodnie z procedurą w Konstrukcji macierzy sztywności stycznej.
  2. Aby nałożyć geometryczne warunki brzegowe, zmodyfikuj macierz sztywności stycznej i wektor residuum dla stopni swobody objętych więzami przemieszczenia, uzyskując \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\) (Obsługa geometrycznych warunków brzegowych).
  3. Rozwiąż równanie liniowe \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\), aby uzyskać poprawkę \(d\boldsymbol{u}_i\). Procedura ta często odpowiada za większość kosztu obliczeniowego iteracji.
  4. Zaktualizuj przyrost przemieszczenia jako \(\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i\), a następnie oblicz odpowiednio wektor sił wewnętrznych \(\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) i residuum \(\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\).
  5. Sprawdź zbieżność i zakończ iterację, jeśli została osiągnięta. Składowe odpowiadające reakcjom więzów pojawiają się w residuum \(\boldsymbol{R}_i\) na stopniach swobody objętych geometrycznymi warunkami brzegowymi, dlatego wskaźnik zbieżności konstruuje się z \(\tilde{\boldsymbol{R}}_i\) po wykluczeniu tych składowych. Konkretne wskaźniki zbieżności i progi opisano w Kryteriach zbieżności. Jeśli zbieżność nie zostanie osiągnięta i osiągnięty zostanie limit iteracji, iterację uznaje się za nieudaną.

Gdy iteracja osiągnie zbieżność, dodaj zbieżne \(\Delta\boldsymbol{u}\) do \(\boldsymbol{u}_n\), aby uzyskać skumulowane przemieszczenie w chwili \(t_{n+1}\), \(\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}\), i przejdź do następnego kroku czasowego.

Tematy powiązane