Aller au contenu

Méthode de Newton-Raphson

Linéarisation et relation de récurrence itérative

Travail virtuel des forces externes et assemblage de l'équation globale donne une équation non linéaire pour le déplacement nodal au temps \(t_{n+1}\), \(\boldsymbol{u}_{n+1}\), qui est résolue par la méthode de Newton-Raphson. Le déplacement nodal jusqu'au temps \(t_n\), \(\boldsymbol{u}_n\), est supposé connu, et l'incrément de déplacement \(\Delta\boldsymbol{u}\) est pris comme variable inconnue à déterminer

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

Dans la suite, la dépendance du vecteur des forces externes vis-à-vis du déplacement nodal est négligée et, avec \(\boldsymbol{F}(\boldsymbol{u}_{n+1}) = \boldsymbol{F}_{n+1}\),

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

est résolue.

Pour la solution actuelle \(\Delta\boldsymbol{u}\), on définit la raideur tangente

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

En l'utilisant, la linéarisation de l'équation non linéaire donne

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

Soit \(d\boldsymbol{u}_i\) la correction à la \(i\)-ème itération, et soit le vecteur résiduel au début de l'itération

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

La relation de récurrence itérative est alors

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

Ainsi, le résidu \(\boldsymbol{R}_i\) est une grandeur correspondant au déséquilibre des forces par rapport à l'équilibre.

Construction de la matrice de raideur tangente

La raideur tangente \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\) est construite en dérivant partiellement par rapport au déplacement nodal le vecteur élémentaire des forces internes obtenu dans Discrétisation du travail virtuel des forces internes, en intégrant les intégrandes obtenues au niveau élémentaire sur chaque domaine d'élément, puis en les assemblant. En notant l'intégrande élémentaire \(\boldsymbol{K}^e_X\) (notation dans la configuration de référence, formulation TL) ou \(\boldsymbol{K}^e_x\) (notation dans la configuration actuelle, formulation UL), la raideur tangente élémentaire est

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

Les formes finales des intégrandes TL/UL sont données ci-dessous. Dans les deux cas, elles sont décomposées en la somme d'un terme de raideur matériau (terme de déplacement initial) et d'un terme de raideur géométrique (terme de contrainte initiale).

Formulation lagrangienne totale

Dans la formulation lagrangienne totale, on suppose une relation linéaire entre le taux de la seconde contrainte de Piola-Kirchhoff \(\dot{\boldsymbol{S}}\) et le taux de déformation de Green-Lagrange \(\dot{\boldsymbol{E}}\), à savoir \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\). Cela correspond aux lois de comportement des matériaux élastiques linéaires (matériaux de St. Venant-Kirchhoff) et des matériaux hyperélastiques, et FrontISTR utilise la formulation lagrangienne totale pour ces matériaux. L'intégrande de la raideur tangente élémentaire s'écrit alors sous forme tensorielle comme

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

Le premier terme du membre de droite est le terme de raideur matériau (terme de déplacement initial), et le second est le terme de raideur géométrique (terme de contrainte initiale).

Dans l'implémentation de FrontISTR, cette intégrande est évaluée sous forme matricielle à l'aide de la notation de 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 \]

Les différentes matrices sont les suivantes. \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}\) sont les matrices B introduites dans Discrétisation du travail virtuel des forces internes, et \(\tilde{\boldsymbol{C}}\) est la représentation de Voigt du tenseur constitutif \(\boldsymbol{\mathsf{C}}\), c'est-à-dire une matrice de raideur matériau \(6\times 6\) (Notation tensorielle et fondements mathématiques). \(\boldsymbol{S}_9, \boldsymbol{F}_9\) sont les matrices de réarrangement suivantes, utilisées pour exprimer le terme de raideur géométrique sous forme d'un produit matriciel. Tout d'abord, pour un tenseur du second ordre \(3\times 3\) \(\boldsymbol{A}\), on définit la notation \([\,\cdot\,]\) qui le réarrange en un vecteur à 9 composantes comme

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

Avec cette définition, \(\boldsymbol{F}_9\) exprime la variation du gradient de déformation sous la forme \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\) et est une matrice \(9\times d n_e\). Pour le nœud d'élément \(\alpha = 1, \ldots, n_e\), le bloc \(9\times d\) correspondant est

\[ [\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{ est la matrice identité } 3\times 3) \]

et \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\), les blocs étant disposés horizontalement dans l'ordre des nœuds de l'élément. \(\boldsymbol{S}_9\) est choisie de telle sorte que, lorsqu'elle est combinée à cette matrice, le terme de raideur géométrique s'exprime sous la forme \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\) ; il s'agit de la matrice \(9\times 9\) suivante

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

C'est la matrice obtenue.

Formulation lagrangienne actualisée

Dans la formulation lagrangienne actualisée, on suppose une relation linéaire entre le taux de Jaumann du tenseur relatif des contraintes de Kirchhoff \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) et le tenseur des taux de déformation \(\boldsymbol{D}\), à savoir \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\). Il s'agit de la forme d'une loi de comportement hypoélastique commune aux matériaux élastiques linéaires, élastoplastiques et de fluage, et FrontISTR utilise la formulation lagrangienne actualisée pour ces matériaux. L'intégrande de la raideur tangente élémentaire exprimée dans la configuration actuelle s'écrit alors sous forme tensorielle comme

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

où \(\boldsymbol{\sigma}^{\nabla T}\) est le taux de Truesdell, \(\boldsymbol{A}_{(L)}\) est la partie linéaire de la déformation d'Almansi, \(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\) est le gradient de déplacement par rapport à la configuration actuelle, et \(\boldsymbol{L}\) est le tenseur gradient de vitesse. Le premier terme du membre de droite est le terme de raideur matériau, et le second est le terme de raideur géométrique.

Dans l'implémentation de FrontISTR, cette intégrande est évaluée sous forme matricielle à l'aide de la notation de 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 \]

Ici, \(\boldsymbol{b}\) est la matrice B construite dans la configuration actuelle (Discrétisation du travail virtuel des forces internes). \(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\) sont obtenues à partir de \(\boldsymbol{S}_9, \boldsymbol{F}_9\) définies pour la formulation TL en remplaçant la seconde contrainte PK \(\boldsymbol{S}\) par la contrainte de Cauchy \(\boldsymbol{\sigma}\) et le gradient dans la configuration de référence \(\partial N_\alpha^e/\partial X_i\) par le gradient dans la configuration actuelle \(\partial N_\alpha^e/\partial x_i\).

\(\boldsymbol{G}\) est une matrice de correction dépendant de la contrainte de Cauchy, nécessaire pour rendre la loi de comportement hypoélastique \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) cohérente avec le cadre de la raideur tangente en tant que loi de comportement basée sur le taux de Truesdell. Elle est obtenue en disposant les composantes du tenseur d'ordre quatre \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) sous forme de Voigt \(6\times 6\) comme

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

C'est la matrice obtenue.

Assemblage de la matrice de raideur globale

La raideur tangente globale \(\boldsymbol{K}\) est obtenue en divisant chaque raideur élémentaire \(\boldsymbol{K}^e\) en blocs \(d\times d\) \(\boldsymbol{K}^e_{\alpha\beta}\) pour chaque paire de nœuds et en utilisant l'ensemble d'assemblage de tenseurs d'ordre deux \(\mathcal{E}^2(i_g, i_h)\) introduit dans Assemblage des grandeurs physiques élément-nœud :

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

Les valeurs obtenues sont disposées dans une matrice dont la ligne est \(i_g\) et la colonne \(i_h\). Dans l'implémentation, l'ensemble \(\mathcal{E}^2\) n'est pas construit explicitement ; les blocs correspondants sont ajoutés directement dans la boucle sur les éléments. La matrice est carrée, de dimension égale au nombre de degrés de liberté par nœud \(\times\) le nombre total de nœuds \(n_g\), mais, les composantes autres que celles reliant des nœuds connectés par des éléments étant \(0\), elle est stockée sous forme de matrice creuse.

Les matrices de raideur élémentaires des formulations TL et UL ont la même forme, à l'exception du changement de configuration de référence (coordonnées nodales et source utilisée pour construire la matrice B) et de la présence ou de l'absence de la matrice \(\boldsymbol{G}\). FrontISTR implémente donc les deux formulations dans un sous-programme commun.

Algorithme d'itération

En résumé, au début de l'itération, posez \(\Delta\boldsymbol{u} = \boldsymbol{0}\) et calculez le résidu initial \(\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)\). Puis, à la \(i\)-ème itération, effectuez la procédure suivante.

  1. Pour le déplacement actuel \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\), calculez la raideur tangente \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\) selon la procédure de Construction de la matrice de raideur tangente.
  2. Pour imposer les conditions aux limites géométriques, modifiez la matrice de raideur tangente et le vecteur résiduel pour les degrés de liberté soumis à des contraintes de déplacement, afin d'obtenir \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\) (Traitement des conditions aux limites géométriques).
  3. Résolvez l'équation linéaire \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\) pour obtenir la correction \(d\boldsymbol{u}_i\). Cette procédure représente souvent la majeure partie du coût de calcul de l'itération.
  4. Mettez à jour l'incrément de déplacement selon \(\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i\), puis calculez en conséquence le vecteur des forces internes \(\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) et le résidu \(\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\).
  5. Vérifiez la convergence et terminez l'itération si elle est atteinte. Les composantes correspondant aux réactions de contrainte apparaissent dans le résidu \(\boldsymbol{R}_i\) aux degrés de liberté soumis à des conditions aux limites géométriques ; l'indicateur de convergence est donc construit à partir de \(\tilde{\boldsymbol{R}}_i\) après exclusion de ces composantes. Les indicateurs et seuils de convergence spécifiques sont décrits dans Critères de convergence. Si la convergence n'est pas atteinte et que la limite d'itérations est atteinte, l'itération est considérée comme ayant échoué.

Lorsque l'itération converge, ajoutez le \(\Delta\boldsymbol{u}\) convergé à \(\boldsymbol{u}_n\) pour obtenir le déplacement cumulé au temps \(t_{n+1}\), \(\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}\), puis passez au pas de temps suivant.

Sujets connexes