Table of Contents

Commit 2016-04-25

Modifications principales

Ce commit porte principalement sur l'amélioration du Lagrangien Augmenté.

Algorithme du Lagrangien Augmenté

Je vous met à profit les derniers développements sur l'algorithme du Lagrangien Augmenté.

Algorithme Général

J'ai constaté à ma grande surprise qu'il faut forcer le calcul d'une itération mécanique après une augmentation. Après chaque augmentation (modification des lagrangiens augmentés à configuration cinématique fixée), on résout l'équilibre mécanique par la méthode de Newton Raphson. Cependant, il se peut que la modification du vecteur des forces externes ne soit pas suffisamment conséquente pour changer la norme moyenne du résidu (alors que la norme max du résidu est nettement plus sensible à ces changements là …) et dès lors forcer une nouvelle itération mécanique, c'est à dire une correction des positions et obtenir une nouvelle configuration cinématique. Donc, l'algorithme continue de modifier les lagrangiens augmentés à configuration cinématique donnée jusqu'à espérer avoir une variation suffisante des forces externes pour avoir une itération mécanique. Si la configuration donnée ne change pas, les gaps normal et tangentiel ne peuvent pas évoluer et nous ne pouvons pas converger vers les tolérances imposées à l'algorithme du lagrangien augmenté. En résumé, l'algorithme stagne et l'intégration temporelle risque de s'arrêter tout simplement !

Il est vrai que les articles de référence en la matière travaillent en général en terme de norme énergétique (ce qui nécessite d'office le calcul d'une itération mécanique !).

Dans le jeux de données Python, on a

alm = AutomaticAugmentedLagrangianManager(metafor)
alm.setForceFirstIterationAfterAugmentation(True)
alm = StandardAugmentedLagrangianManager(metafor)
alm.setForceFirstIterationAfterAugmentation(True)

Par défaut, ce paramètre est mis à False !

Mise à jour des lagrangiens augmentés

J'ai développé une autre approche pour mettre à jour les lagrangiens augmentés sur base de l'histoire de l'évolution de la pression de contact en fonction du gap normal et le cisaillement de contact en fonction du gap tangentiel. Ces lagrangiens augmentés “extrapolés” sont mis à jour si besoin via le critère de frottement de Coulomb ou de Tresca. Au lieu de mettre à jour les coefficients de pénalité pour chaque nœud de contact à chaque augmentation ou au cours de chaque itération mécanique, on utilise cette autre approche développée par Wriggers and Zavarise en contact sans frottement qui possède un ordre de convergence super-linéaire, mais qui n'est plus spécialement monotone en convergence sans garde fou !

Avec ces derniers gardes fous, l'algorithme semble être encore plus efficace que l'approche standard car on tient de l'historique au cours des augmentations !

augLagAugmentation = AugLagExtrapolationAugmentation(alm)
augLagAugmentation.setUseStabilityCriterion(True)
augLagAugmentation.setUseExtrapolationCriterion(True)
augLagAugmentation.setExtrapolationFactor(_factor)
augLagAugmentation.setResetPreviousAugmentedLagrangians(True)

La dernière option permet d'activer la procédure dès la première augmentation en indiquant False, sinon il faut au moins une augmentation selon la procédure d'Uzawa avant de lancer celle par extrapolation.

Par défaut, tous les critères sont activés et le facteur d'extrapolation est 0.1 (c'est à dire on autorise une variation de maximum 10 % de la valeur courante de la pression de contact et du cisaillement de contact.).

Cette procédure est activée si on a le nœud de contact est en contact au cours de deux résolutions successives de l'équilibre mécanique (Le lagrangien normal est mis à jour en fonction de la valeur à l'augmentation k du couple (contact pressure, normal gap) et l'augmentation k-1 du même couple), sinon les lagrangiens augmentés sont mis à jour selon la méthode d'Uzawa.

Finalement, il est possible de mettre à zéro les lagrangiens augmentés d'un pas de temps à l'autre (très déconseillés !), quelque soit la méthode d'augmentation :

augLagAugmentation.setResetAugmentedLagrangians(True)

Critère d'arrêt géométrique

J'ai ajouté une nouvelle méthode de calcul de longueur de contact utilisée lors des critères géométriques normalisés. Dans les estimateurs de pas de temps critique dans les schémas explicites, la notion de dimension équivalente d'un élément apparaît naturellement dans la définition du pas de temps critique. Cependant, ce choix reste vague (distance minimale entre deux nœuds de la maille ou autre) et dès lors, il propose une analyse des valeurs propres à partir de laquelle il dérive une longueur équivalente (plus conservative !). Leur raisonnement est indépendant à priori du type de maille et de la dimension de l'espace. C'est pour ça, que j'ai repris leur développement et je l'ai codé dans Metafor pour le cas 2D, 3D et axisymmétrique.

En résumé, on a les méthodes suivantes :

Par exemple, pour un élément quadrangulaire rectangulaire $a \times b$ (avec $b > a > 0$), on a

Dans le jeux de données Python, on choisit les options comme ceci :

augLagGeoCriterion = AugLagNormalisedGeoCriterion(alm)
augLagGeoCriterion.setCharacteristicLengthMeth(_method)

Sur un pas de temps donné, la maille peut subir de grande déformation d'où le fait que l'on recalcule à chaque fois la longueur caractéristique au cours des augmentations.

Mise à jour du pas de temps

L'algorithme du Lagrangien augmenté permet à priori de résoudre un problème de contact sur un pas de temps plus élevé que l'équivalent en Méthode de Pénalité. Pour mettre à jour le pas de temps, j'ai proposé plusieurs méthodes :

Cette valeur d'itération “équivalente” est insérée dans l'algorithme de mise à jour du pas de temps proposé par JPP dans sa thèse.

Dans le jeux de données Python, on a

tscm = NbOfAugLagMechNRIterationsTimeStepComputationMethod(metafor)
tscm.setTimeStepDivisionFactor(_timeStepDivisionFactor)
tscm.setEquivalentIteMethod(_equivalentIteMethod)
tscm.setMaxNbOfIteInAdaptOfTimeStep(_maxNbOfIteInAdaptOfTimeStep)
tscm.setNbOptiIte(_nbOptiIte)
tsm.setTimeStepComputationMethod(tscm)

Par exemple, si on a le profil suivant de convergence sur un pas de temps de donné (12 augmentations pour 30 itérations mécaniques) :

ALM ITE NO MECH ITE NUMBER
0 2
1 3
2 4
3 3
4 3
5 3
6 2
7 2
8 2
9 2
10 2
11 1
12 1

on obtient l'itération mécanique équivalente suivante :

ALMTSCM EQUIVALENT MECH ITE NUMBER
ALM_TSCM_FIRST 2
ALM_TSCM_MAX 4
ALM_TSCM_MEAN 2
ALM_TSCM_WEIGHTEDMEAN 3

Informations sur Interaction de Contact

Il est possible de demander d'afficher des informations sur le status des nœuds de contact d'une interaction.

ci = RdContactInteraction(1)
ci.setShowInformation(True)

Voici un exemple de sortie.

mechanical iteration 0  : Rmean = 1.000000e+000  Rmax = 2.244123e+002 
RdContactInteraction #2
Actual/Total Number of Nodes in Contact = 4 / 11
Number of Nodes in Sticking/Sliding Contact = 4 / 0
mechanical iteration 1  : Rmean = 2.192696e-001  Rmax = 4.097902e+001 
RdContactInteraction #2
Actual/Total Number of Nodes in Contact = 4 / 11
Number of Nodes in Sticking/Sliding Contact = 4 / 0
mechanical iteration 2  : Rmean = 2.016662e-001  Rmax = 2.821996e+001 
RdContactInteraction #2
Actual/Total Number of Nodes in Contact = 3 / 11
Number of Nodes in Sticking/Sliding Contact = 3 / 0
mechanical iteration 3  : Rmean = 1.084675e-001  Rmax = 9.532517e+000 
RdContactInteraction #2
Actual/Total Number of Nodes in Contact = 3 / 11
Number of Nodes in Sticking/Sliding Contact = 3 / 0
mechanical iteration 4  : Rmean = 2.831632e-002  Rmax = 2.978541e+000 
RdContactInteraction #2
Actual/Total Number of Nodes in Contact = 3 / 11
Number of Nodes in Sticking/Sliding Contact = 3 / 0
mechanical iteration 5  : Rmean = 1.060746e-003  Rmax = 1.218047e-001 
RdContactInteraction #2
Actual/Total Number of Nodes in Contact = 3 / 11
Number of Nodes in Sticking/Sliding Contact = 3 / 0
mechanical iteration 6  : Rmean = 1.744038e-005  Rmax = 1.792596e-003 
RdContactInteraction #2
Actual/Total Number of Nodes in Contact = 3 / 11
Number of Nodes in Sticking/Sliding Contact = 3 / 0

Cela peut permettre de comprendre la raison pour laquelle un pas de temps ne converge pas (Status collant glissant qui change ou status en contact ou hors contact qui change).

Nombre de conditionnement réciproque

Il est désormais possible d'estimer le nombre de conditionnement réciproque d'une matrice carrée (Matrix de Metafor) quelconque via les routines LAPACK. Les matrices utilisées pour la résolution du système linéaire sous-jacente au calcul de l'équilibre implicite se trouvent sous un forme creux liée au solveur utilisé. Cependant, il n'existe pas (sauf erreur de ma part) d'estimateur du nombre de conditionnement réciproque dans les routines LAPACK, donc nous devons transformer une matrice creuse sous un format plein pour pouvoir effectuer le calcul proprement dit (ce n'est pas top du tout mais bon faute de mieux …). J'ai réuni un peu de documentation sur ces estimateurs, si qqn veut avoir les infos.

On peut demander de calculer à chaque résolution du système linéaire, le nombre de conditionnement réciproque. Pour se faire, il suffit de faire :

solvermanager = metafor.getSolverManager()
solver = DSSolver()
solver.setShowRCond(True)
solvermanager.setSolver(solver)        

Également, il est possible de faire un extracteur de cette grandeur (Nombre de conditionnement réciproque de la matrice FreeJacobian du schéma d'intégration temporelle) via les commandes suivantes :

valuesmanager.add(1, MiscValueExtractor(metafor, EXT_RCOND), 'RCond')

Dans Matlab, on compare ce nombre de conditionnement réciproque avec l'epsilon machine multipliée par la taille du système à résoudre. Pour rappel, ce nombre est proche de 1 lorsque le système est bien conditionné et ce nombre est proche de 0 lorsque le système est mal conditionné.

Finalement, cet extracteur peut aider à quantifier l'influence du coefficient de pénalité sur la résolution du système linéaire.

Nouveaux extracteurs de contact

Par la même occasion, j'ai ajouté plein de nouveau extracteurs de variables de contact. J'ai ajouté directement les infos dans la partie sur les extracteurs de contact. Finalement, la documentation est à jour à ce niveau là.

Fichiers/Dossiers ajoutés/supprimés

[a]:mtFEMBase/AugLagTimeStepComputationMethod.cpp
[a]:mtFEMBase/AugLagTimeStepComputationMethod.h
[a]:mtContact/src/AugLagExtrapolationAugmentation.cpp
[a]:mtContact/src/AugLagExtrapolationAugmentation.h
[a]:mtContact/src/AugLagNormalExtrapolationAugmentation.cpp
[a]:mtContact/src/AugLagNormalExtrapolationAugmentation.h
[a]:mtContact/src/AugLagTangentialExtrapolationAugmentation.cpp
[a]:mtContact/src/AugLagTangentialExtrapolationAugmentation.h
[a]:mtContact/src/AugLagFrictionalExtrapolationAugmentation.cpp
[a]:mtContact/src/AugLagFrictionalExtrapolationAugmentation.h
[a]:mtContact/src/NormalAugmentedLagrangianValueExtractor.cpp
[a]:mtContact/src/NormalAugmentedLagrangianValueExtractor.h
[r]:mtContact/src/GapValueExtractor.cpp => [a]:mtContact/src/NormalGapValueExtractor.cpp
[r]:mtContact/src/GapValueExtractor.h => [a]:mtContact/src/NormalGapValueExtractor.h
[a]:mtContact/src/Shear1ContactValueExtractor.cpp
[a]:mtContact/src/Shear1ContactValueExtractor.h
[a]:mtContact/src/Shear2ContactValueExtractor.cpp
[a]:mtContact/src/Shear2ContactValueExtractor.h
[a]:mtContact/src/Tangent1ForceValueExtractor.cpp
[a]:mtContact/src/Tangent1ForceValueExtractor.h
[a]:mtContact/src/Tangent2ForceValueExtractor.cpp
[a]:mtContact/src/Tangent2ForceValueExtractor.h
[a]:mtContact/src/TangentAugmentedLagrangian1ValueExtractor.cpp
[a]:mtContact/src/TangentAugmentedLagrangian1ValueExtractor.h
[a]:mtContact/src/TangentAugmentedLagrangian2ValueExtractor.cpp
[a]:mtContact/src/TangentAugmentedLagrangian2ValueExtractor.h
[a]:mtContact/src/TangentAugmentedLagrangianValueExtractor.cpp
[a]:mtContact/src/TangentAugmentedLagrangianValueExtractor.h
[a]:mtContact/src/TangentGap1ValueExtractor.cpp
[a]:mtContact/src/TangentGap1ValueExtractor.h
[a]:mtContact/src/TangentGap2ValueExtractor.cpp
[a]:mtContact/src/TangentGap2ValueExtractor.h
[a]:mtContact/src/TangentGapValueExtractor.cpp
[a]:mtContact/src/TangentGapValueExtractor.h

[r]:

Tests ajoutés/supprimés

[r]:
[a]:apps/qs/cont2AugLag.py

Gaëtan WAUTELET 2016/04/25