Vai al contenuto

Metodo di Newton-Raphson

Linearizzazione e relazione di ricorrenza iterativa

Il lavoro virtuale delle forze esterne e l'assemblaggio dell'equazione globale fornisce un'equazione non lineare per lo spostamento nodale al tempo \(t_{n+1}\), \(\boldsymbol{u}_{n+1}\), che viene risolta mediante il metodo di Newton-Raphson. Si assume noto lo spostamento nodale fino al tempo \(t_n\), \(\boldsymbol{u}_n\), e si considera l'incremento di spostamento \(\Delta\boldsymbol{u}\) come variabile incognita da determinare

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

Nel seguito si trascura la dipendenza del vettore delle forze esterne dallo spostamento nodale e, ponendo \(\boldsymbol{F}(\boldsymbol{u}_{n+1}) = \boldsymbol{F}_{n+1}\),

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

si risolve tale equazione.

Alla soluzione corrente \(\Delta\boldsymbol{u}\) si definisce la rigidezza tangente

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

Usandola, la linearizzazione dell'equazione non lineare fornisce

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

Sia \(d\boldsymbol{u}_i\) la correzione alla \(i\)-esima iterazione e sia il vettore residuo all'inizio dell'iterazione

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

La relazione di ricorrenza iterativa è quindi

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

Pertanto il residuo \(\boldsymbol{R}_i\) è una quantità che corrisponde allo squilibrio delle forze rispetto alla condizione di equilibrio.

Costruzione della matrice di rigidezza tangente

La rigidezza tangente \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\) viene costruita derivando parzialmente rispetto allo spostamento nodale il vettore delle forze interne dell'elemento ottenuto in Discretizzazione del lavoro virtuale delle forze interne, integrando su ciascun dominio di elemento gli integrandi risultanti a livello di elemento e assemblandoli. Indicando l'integrando a livello di elemento con \(\boldsymbol{K}^e_X\) (notazione nella configurazione di riferimento, formulazione TL) oppure \(\boldsymbol{K}^e_x\) (notazione nella configurazione corrente, formulazione UL), la rigidezza tangente dell'elemento è

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

Di seguito sono riportate le forme finali degli integrandi TL/UL. In entrambi i casi essi sono decomposti nella somma di un termine di rigidezza del materiale (termine di spostamento iniziale) e di un termine di rigidezza geometrica (termine di tensione iniziale).

Formulazione Total Lagrange

Nella formulazione Total Lagrange si assume una relazione lineare tra la velocità della seconda tensione di Piola-Kirchhoff \(\dot{\boldsymbol{S}}\) e la velocità della deformazione di Green-Lagrange \(\dot{\boldsymbol{E}}\), ossia \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\). Ciò corrisponde alle leggi costitutive dei materiali elastici lineari (materiali di St. Venant-Kirchhoff) e dei materiali iperelastici; FrontISTR utilizza la formulazione Total Lagrange per tali materiali. L'integrando della rigidezza tangente dell'elemento viene quindi scritto in forma tensoriale come

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

Il primo termine a destra è il termine di rigidezza del materiale (termine di spostamento iniziale), mentre il secondo è il termine di rigidezza geometrica (termine di tensione iniziale).

Nell'implementazione di FrontISTR, questo integrando viene valutato in forma matriciale mediante la notazione di Voigt:

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

Le singole matrici sono le seguenti. \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}\) sono le matrici B introdotte in Discretizzazione del lavoro virtuale delle forze interne, mentre \(\tilde{\boldsymbol{C}}\) è la rappresentazione di Voigt del tensore costitutivo \(\boldsymbol{\mathsf{C}}\), cioè una matrice di rigidezza del materiale \(6\times 6\) (Notazione tensoriale e fondamenti matematici). \(\boldsymbol{S}_9, \boldsymbol{F}_9\) sono le seguenti matrici di riordinamento usate per esprimere il termine di rigidezza geometrica come prodotto matriciale. Innanzitutto, per un tensore del secondo ordine \(3\times 3\) \(\boldsymbol{A}\), si definisce la notazione \([\,\cdot\,]\) che lo riordina in un vettore a 9 componenti come

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

Con questa definizione, \(\boldsymbol{F}_9\) esprime la variazione del gradiente di deformazione nella forma \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\) ed è una matrice \(9\times d n_e\). Per il nodo di elemento \(\alpha = 1, \ldots, n_e\), il corrispondente blocco \(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{ è la } 3\times 3 \text{ matrice identità}) \]

ed è dato da \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\), con i blocchi disposti orizzontalmente nell'ordine dei nodi dell'elemento. \(\boldsymbol{S}_9\) è scelta in modo che, in combinazione con questa matrice, il termine di rigidezza geometrica sia espresso come \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\); è la seguente matrice \(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} \]

Questa è la matrice risultante.

Formulazione Updated Lagrange

Nella formulazione Updated Lagrange si assume una relazione lineare tra la velocità di Jaumann del tensore di tensione relativa di Kirchhoff \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) e il tensore velocità di deformazione \(\boldsymbol{D}\), ossia \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\). Questa è la forma di una legge costitutiva ipoelastica comune ai materiali elastici lineari, elastoplastici e con creep; FrontISTR utilizza la formulazione Updated Lagrange per tali materiali. L'integrando della rigidezza tangente dell'elemento espresso nella configurazione corrente viene quindi scritto in forma tensoriale come

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

dove \(\boldsymbol{\sigma}^{\nabla T}\) è la velocità di Truesdell, \(\boldsymbol{A}_{(L)}\) è la parte lineare della deformazione di Almansi, \(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\) è il gradiente di spostamento rispetto alla configurazione corrente e \(\boldsymbol{L}\) è il tensore gradiente di velocità. Il primo termine a destra è il termine di rigidezza del materiale, mentre il secondo è il termine di rigidezza geometrica.

Nell'implementazione di FrontISTR, questo integrando viene valutato in forma matriciale mediante la notazione di Voigt:

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

Qui \(\boldsymbol{b}\) è la matrice B costruita nella configurazione corrente (Discretizzazione del lavoro virtuale delle forze interne). \(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\) si ottengono da \(\boldsymbol{S}_9, \boldsymbol{F}_9\) definite per la formulazione TL sostituendo la seconda tensione PK \(\boldsymbol{S}\) con la tensione di Cauchy \(\boldsymbol{\sigma}\) e il gradiente nella configurazione di riferimento \(\partial N_\alpha^e/\partial X_i\) con il gradiente nella configurazione corrente \(\partial N_\alpha^e/\partial x_i\).

\(\boldsymbol{G}\) è una matrice di correzione dipendente dalla tensione di Cauchy, necessaria per rendere la legge costitutiva ipoelastica \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) coerente con il quadro della rigidezza tangente formulato mediante una legge costitutiva basata sulla velocità di Truesdell. Si ottiene disponendo le componenti del tensore del quarto ordine \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) nella forma di Voigt \(6\times 6\) come

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

Questa è la matrice risultante.

Assemblaggio della matrice di rigidezza globale

La rigidezza tangente globale \(\boldsymbol{K}\) si ottiene suddividendo ciascuna rigidezza di elemento \(\boldsymbol{K}^e\) in blocchi \(d\times d\) \(\boldsymbol{K}^e_{\alpha\beta}\) per ciascuna coppia di nodi e utilizzando l'insieme di assemblaggio per tensori del secondo ordine \(\mathcal{E}^2(i_g, i_h)\) introdotto in Assemblaggio delle grandezze fisiche nodali degli elementi:

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

I valori risultanti sono disposti in una matrice con riga \(i_g\) e colonna \(i_h\). Nell'implementazione, l'insieme \(\mathcal{E}^2\) non viene costruito esplicitamente; i blocchi corrispondenti vengono invece sommati direttamente all'interno del ciclo sugli elementi. La matrice è quadrata e ha dimensione pari ai gradi di libertà per nodo \(\times\) il numero totale di nodi \(n_g\), ma poiché le componenti diverse da quelle tra nodi collegati attraverso gli elementi sono pari a \(0\), viene memorizzata in forma sparsa.

Le matrici di rigidezza degli elementi per le formulazioni TL e UL hanno la stessa forma, salvo la commutazione della configurazione di riferimento (coordinate nodali e sorgente usata per costruire la matrice B) e la presenza o assenza della matrice \(\boldsymbol{G}\). FrontISTR implementa quindi entrambe le formulazioni in una subroutine comune.

Algoritmo iterativo

Riassumendo quanto sopra, all'inizio dell'iterazione si imposta \(\Delta\boldsymbol{u} = \boldsymbol{0}\) e si calcola il residuo iniziale \(\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)\). Quindi, alla \(i\)-esima iterazione, si esegue la procedura seguente.

  1. Allo spostamento corrente \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\), calcolare la rigidezza tangente \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\) mediante la procedura descritta in Costruzione della matrice di rigidezza tangente.
  2. Per imporre le condizioni al contorno geometriche, modificare la matrice di rigidezza tangente e il vettore residuo per i gradi di libertà soggetti a vincoli di spostamento, ottenendo \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\) (Trattamento delle condizioni al contorno geometriche).
  3. Risolvere l'equazione lineare \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\) per ottenere la correzione \(d\boldsymbol{u}_i\). Questa procedura rappresenta spesso la maggior parte del costo computazionale del calcolo iterativo.
  4. Aggiornare l'incremento di spostamento come \(\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i\) e, di conseguenza, calcolare il vettore delle forze interne \(\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) e il residuo \(\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\).
  5. Verificare la convergenza e terminare l'iterazione se è stata raggiunta. Le componenti corrispondenti alle reazioni vincolari compaiono nel residuo \(\boldsymbol{R}_i\) nei gradi di libertà soggetti a condizioni al contorno geometriche; pertanto l'indicatore di convergenza viene costruito a partire da \(\tilde{\boldsymbol{R}}_i\) dopo aver escluso tali componenti. Gli indicatori e le soglie di convergenza specifici sono descritti in Criteri di convergenza. Se la convergenza non viene raggiunta e si raggiunge il limite di iterazione, l'iterazione viene considerata fallita.

Quando l'iterazione converge, si somma il \(\Delta\boldsymbol{u}\) convergente a \(\boldsymbol{u}_n\) per ottenere lo spostamento accumulato al tempo \(t_{n+1}\), \(\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}\), e si procede al passo temporale successivo.

Argomenti correlati