Le formalisme Arbitraire Eulérien Lagrangien (ALE) est développé depuis de nombreuses années dans Metafor. Il permet de découpler le déplacement du maillage et le déplacement de la matière.

Les buts recherchés peuvent être:

• minimiser la distorsion du maillage: on sait que la solution est d'autant plus précise que les mailles sont proches de quadrangles/hexaèdres parfaits. Le formalisme ALE permet d'améliorer l'aspect des mailles sans pour autant remailler complètement la structure (cette dernière opération serait non seulement coûteuse en temps de calcul mais aussi diffusive et donc peu précise).
• minimiser la taille du maillage en raffinant les zones qui nécessitent de nombreuses mailles sans se soucier du déplacement de la matière (par exemple garder un grand nombre de mailles près des contacts).
• résoudre des problèmes eulériens ou quasi-eulériens (problèmes stationnaires).

Dans Metafor, on se limite à une résolution découplée du problème (on parle d' “operator split”), ceci veut dire que, pour chaque pas de temps, on résout le problème comme s'il était lagrangien. A la fin du pas de temps, lorsqu'on a convergé, on effectue successivement les opérations suivantes:

• repositionnement des noeuds (rezoning) : metafor dispose d'une série de méthodes pour définir le plus simplement possible le mouvement des noeuds par rapport à celui de la matière.
• convection de grandeurs utiles d'un maillage vers l'autre : actuellement l'algorithme utilisé est un algorithme de convection par volumes finis.
• post-traitement: recalcul des forces après convection, ajustement des contacts, etc

Utiliser le formalisme ALE revient donc à définir des opérateurs de repositionnement de noeuds qui seront appelés à chaque activation de l'algorithme et paramétrer un schéma de convection.

L'objet gérant l'ALE se trouve dans l'analyse Metafor :

ale = AleMethod(metafor)

Le passé nous a montré que l'utilisation de l'ALE est réservé aux experts si un effort d'interfaçage important n'est pas réalisé (cfr. thèse de JPP). Un travail très important a été apporté dans la version actuelle de Metafor pour faciliter son utilisation. Cela reste cependant délicat.

Activer l'ALE se fait donc en 1 commande:

ale.enable(n)

où n est le nombre de pas entre 2 activations de l'algorithme (par défaut n=1). Inversément, supprimer cette ligne d'un test désactivera l'ALE et fera tourner le test en lagrangien.

L'objet ALE contient deux objets représentant la phase de relocalisation des noeuds (ReZoning) et la phase de convection:

rez  = ale.getReZoningStep()
conv = ale.getTransferStep()

J'utilise le terme “rezoning” (repositionnement) et pas “smoothing” (lissage) car le but n'est pas nécessairement un lissage du maillage.

Les commandes de rezoning permettent de définir le déplacement du maillage au cours du temps (éventuellement en fonction de celui de la matière). Chaque méthode est gérée par un objet. Par exemple EulerianReZoner gère un déplacement eulérien du maillage.

Chaque méthode s'applique sur une entité géométrique target (ceci n'est pas restrictif puisqu'un groupe de noeuds est une entité géométrique - la structure ne doit donc pas nécessairement posséder de géométrie).

Ajouter une méthode revient à ajouter un objet de type ReZoner au ReZoningStep de AleMethod.

rez.add(EulerianReZoner(curveset(3)))

Une autre possibilité pour piloter le déplacement du maillage est d'utiliser des “fixations maillage”.

Je liste ici les différentes méthodes de rezoning disponible dans Metafor. Vu la large gamme de problèmes pouvant être traités par Metafor, nous avons besoin d'une large gamme de méthodes de rezoning. Chaque méthode possède des qualités et des défauts: il n'existe pas de méthode universelle. La meilleure méthode dépend souvent de la physique du problème et du type de maillage utilisé (structuré ou non, taille de mailles constantes ou non, etc).

Remarques:

• Les ReZoners sont appliqués dans l'ordre dans lequel ils sont ajoutés au ReZoningStep. Ceci peut avoir de l'influence sur le résultat obtenu.
• Il est très important que l'utilisateur utilise des méthodes qui respectent la physique du problème. En particulier, un noeud sur une surface libre d'un objet doit rester sur cette surface à moins que cette surface ne délimite qu'artificiellement le problème et que de la matière n'entre ou sorte réellement par celle-ci. Aucune vérification de la conservation des surfaces externes d'un solide n'est réalisée! Il est donc possible de déplacer le maillage sans respecter les frontières du problème. Dans ce cas, les résultats n'auront aucun sens!
• Les méthodes surfaciques peuvent être appliquées à des surfaces courbes. Dans ce cas, il faut spécifier si la courbure doit être préservée ou pas (fct membre ReZoner::keepShape()):

Exemple:

rezoner = Tm2DReZoniner(sidset(2))
rezoner.keepShape()
rez.add(rezoner)

Pour l'instant la surface est conservée en projetant sur une surface spline construite sur les patchs arrivant au noeud.

• Pour les méthodes itératives, il est possible d'effectuer une relaxation pour éviter un déplacement brusque du maillage par la commande IterativeReZoner::setRelaxation(alpha).
• Les méthodes itératives prennent un deuxième argument lors de leur construction: le nombre d'itérations à effectuer.

Cette méthode (développée par Tian Zhou) essaye d'égaliser en chaque point les angles que forment les arêtes arrivant à ce point. C'est une méthode 2D.

D'après mes essais, la méthode fonctionne plutôt bien pour des maillages non structurés pour lesquels les mailles ont environ la même taille.

On peut observer des retournements de mailles. Pour un maillage structuré, je conseille d'utiliser cette méthode pondérée avec une autre.

Les noeuds de l'entité cible restent fixes au cours du temps.

Cette méthode est le lissage laplacien (appelée parfois aussi “méthode du barycentre”, “lissage elliptique”, etc ). Les noeuds sont déplacés itérativement au barycentre des centres de gravité des facettes y aboutissant.

Malgré sa simplicité, la méthode est efficace sur des maillages dont la taille des mailles est plus ou moins constante (si ce n'est pas le cas, la méthode va uniformiser les volumes des mailles). On observe cependant des retournements de mailles près des frontières trop convexes.

Plusieurs choix sont possibles:

• Laplacian2DReZoner::useEdges() : utilise les noeuds arêtes (par défaut, utilise le centre de gravité des faces).
• Laplacian2DReZoner::useWeightedForm() : utilise une forme pondérée (par la taille de l'entité choisie - faces ou arête). Par exemple, la combinaison de useEdges() et useWeightedForm() donne la méthode appelée “length weighted laplacian smoothing” et “weighted area method” si on utilise les faces.

Permet de définir une frontière “eulérienne” qui ne sera jamais traversée par le maillage. ceci est très utile pour délimiter un maillage dans l'espace pour la simulation de process stationnaires.

Le contructeur de ce ReZoner prend un argument supplémentaire : la courbe/surface frontière.

Expl:

rezo = BoundaryReZoner(curset(1), srfset(101))

(Vieille) méthode équivalente à la précédente utilisable uniquement dans des tests bidimensionnels pur piloter les points extrêmes d'une ligne frontière. Je conseille d'utiliser BoundaryReZoner à la place.

Optimisation itérative des éléments adjacents à chaque noeud par la méthode de Giuliani. Cette méthode fonctionne mieux que le lissage laplacien mais est plus coûteuse. Elle tend à homogénéiser les volumes des mailles et peut parfois provoquer des retournements de mailles.

Méthode de JPP, utilisable pour remailler des courbes dans le cas de problèmes bi-dimensionnels. Le principe est simple: on fait passer un arc de cercle par le point considéré et les 2 points les plus proches de celui-ci. On replace ensuite le point sur cet arc en respectant l'abscisse curviligne voulue. Méthode 2D uniquement.

Utilise une spline cubique pour remailler les lignes. Fonctionne aussi bien en 2D qu'en 3D. Permet également de remailler des contours (Wires).

Remaille une face ou un volume en utilisant le mailleur transfini. L'objet doit avoir été préalablement maillé avec ce même mailleur.

Appliqué à une surface 3D courbe, il faut spécifier si la courbure doit être préservée ou pas (fct membre ReZoner::keepShape()):

Bidouille de JPP utile pour remailler une ligne droite uniformément entre ses deux points extrêmes.

Cette méthode est une extension du lissage laplacien. On l'appelle aussi lissage elliptique. Elle est applicable pour des maillages structurés de quadrangles. Deux méthodes ont été implémentées : la méthode de Winslow (Equipotential::useWinslow()) et la méthode isoparamétrique (Equipotential::useIsoParametric()).

Ces méthodes sont efficaces sur des maillages structurés lorsque les tailles de mailles varient (elles conservent leur aspect).

Extension du cas précédent (Equipotential::useIsoParametric()) aux maillages surfaciques de quads non structurés. La méthode peut être combinée avec un lissage laplacien par la fonction setWeight(weight).

Cette méthode combine linéairement deux ReZoners avec un facteur de pondération (linear blending).

Expl:

rez1 = LaplacianReZoner(sideset(1), 1)
rez2 = AreaPullReZoner(sideset(1), 1)
rezo = BlendedReZoner(sideset(1), rez1, rez2, 0.9, 10)
ale.getReZoningStep().add(rezo)

combine un “area-pull” et un lissage laplacien classique avec un facteur 0.9 et 10 itérations.

(1) Quads/Hexa only!
(2) Structuré

Il est possible de piloter le maillage à l'aide de fonction définies par un ensemble de points exactement comme c'est le cas pour les fixations “matière”.

Les fixations du maillage sont une manière simple de piloter la position le maillage en ALE en fonction du temps. Il est conseillé d'utiliser cette méthode de pilotage du maillage uniquement lorsque les méthodes de repositionnement ne sont pas satisfaisantes pour le problème traité.

Voir Prescribed Displacements. Il suffit de remplacer le Field classique par un Field maillage ((TX,RRE) = déplacement maillage selon x).

Les fixations "matière" sont des déplacements imposés sur la matière. Dans le cas de l'ALE, ces fixations sont appliquées sur le maillage pendant la phase lagrangienne de l'ALE (voir “operator split”). Lorsqu'une ligne n'est pas lagrangienne, il est souvent utile de pouvoir appliquer ces déplacements de manière incrémentale. Cela veut dire qu'à un pas de temps donné, on déplace le maillage de l'incrément de déplacement. Pour obtenir un calcul incrémental on utilise l'option INCREMENTAL_LOAD lors de la création du déplacement. Par exemple:

domain.getLoadingSet().define(curset(1), Field1D(TX,RRE), 1.0, fct, INCREMENTAL_LOAD )

applique un déplacement incrémental suivant la fontion fct selon x (TX,RRE) sur les noeuds de la courbe 1.

Note: depuis 2008, INCREMENTAL_LOAD est devenu le défaut et peut donc être supprimé de la commande.

La deuxième étape de la deuxième phase de l'algorithme ALE est la convection des données d'un maillage à l'autre.

De quelles données s'agit-il? Ces données sont fonctions de plusieurs facteurs. Elles dépendent:

• du matériau : il faut par exemple transférer les contraintes si le matériau est mécanique, les défos plastiques si le matériau est élastoplastique, etc
• du type d'élément : on doit transférer la pression séparément des grandeurs déviatoriques lorsqu'on utilise un élément sous-intégré (type SRI).
• du schéma d'intégration : si le schéma est dynamique, il faut transférer la densité, les vitesses et éventuellement les accélérations si le schéma est implicite.

Elles sont donc de deux types: des données “aux points d'intégration” (comme par exemple les contraintes) et les données “nodales”.

Tous ces choix sont automatiques et l'utilisateur ne doit pas s'en préoccuper.

Pour rappel, l'objet gérant la convection est accessible via:

conv = ale.getTransferStep()

Pour pouvoir faire des problèmes mixtes où seule une partie du problème est ALE, on définit des “régions”. Chaque région peut être configurée séparement. En particulier, on peut définir des convecteurs différents sur chaque région.

Pour simplifier les choses, chaque interaction de Metafor peut engendrer une région ALE. Il n'est donc pas possible de faire 2 régions hors d'une interaction. Ceci n'est pas restrictif puisqu'une interaction peut facilement être découpée en 2. L'avantage de cette manière de faire est qu'il est alors possible de faire l'hypothèse d'un seul matériau et un seul type d'élément par région.

Pour définir une région ALE à partir d'une interaction noméée app, on écrit:

region = conv.get(app)

La région va générer les convecteurs nécessaires en fonction des éléments, du matériau et du schéma d'intégration. A ce niveau, on peut spécifier des champs que la région doit ignorer lors de la génération des convecteurs par la commande ignore:

region.ignore(IF_FTOTAL)

Dans ce cas, on veut ignorer la convection du tenseur F. ceci aura pour conséquence que toutes les mesures de défos seront fausses en sortie. Si on n'a pas besoin de ces grandeurs, on économise 9 convections! Cette commande est dangereuse puisqu'il est ainsi possible de désactiver par exemple le transfert des contraintes (ce qui, dans tous les cas, entrainera une solution erronée).

Il existe actuellement deux types de convecteurs. Tout deux utilisent la méthode des volumes finis et créent des maillages auxiliaires pour gérer le transfert des données. Ces maillages sont générés automatiquement en fonction des grandeurs à convecter. A chaque maillage est associé un convecteur et une liste de variables.

Le type de convection est géré par le type de l'élément de convection utilisé:

• GodunovConvCell : schéma de Godunov (premier ordre sur maillage structuré)
• LinearRecConvCell : schéma VF à reconstruction linéaire (second ordre sur maillage structuré)

Le choix de la méthode est fait par ajout d'un ElementProperties à la région ALE. Par exemple:

cPrp = ElementProperties(GodunovConvCell)
cPrp.put (UPWIND_COEFF, 1.0)
region.add(cPrp)

De cette manière, tous les convecteurs générés sur cette région utiliseront ces propriétés. Il est également possible de configurer chaque convecteur séparement. pour les désigner, il suffit de spécifier une variable:

region.add(IF_DEV_SIG_XX, cPrp)

Le convecteur qui transfère les contraintes et toutes les autres grandeurs définies sur les mêmes points d'intégration que les contraintes utilisera les propriétés de cPrp.

L'utilisateur peut choisir un coefficient d'upwind (UPWIND_COEFF). Celui-ci sera utilisé pour calculer le flux entre les cellules. Une valeur de 1.0 (par défaut) donne un schéma appellé “full donor cell differencing” (correspond à un décentrage arrière complet). Il est fort diffusif mais est TVD (total variation diminishing - pas d'oscillations) et très rapide.

Cette méthode plus précise mais plus coûteuse permet de capter de manière plus correcte d'éventuelle discontinuités (gradients élevés) et montre beaucoup moins de diffusion. Elle n'est cependant pas TVD si on n'utilise pas un limiteur de flux.

L'utilisateur peut choisir un type de stencil à utiliser (STENCILTYPE) et un limiteur de flux (LIMITERTYPE).

• Le type de stencil:
  • LEASTSQUARE_STENCIL : calcul du gradient par moindre carrés
  • GREENGAUSS_STENCIL : calcul du gradient par la formule de Green-Gauss
• Le type de limiteur:
  • BARTH_LIMITER : limite le flux par l'algorithme de Barth
  • ANISOTROPIC_LIMITER : limite le flux par l'algorithme de Barth modifié de manière anisotrope

Exemple:

cPrp = ElementProperties(LinearRecConvCell)
cPrp.put (STENCILTYPE, LEASTSQUARE_STENCIL)
cPrp.put (LIMITERTYPE, BARTH_LIMITER)
region.add(cPrp)

Les méthodes de convection étant explicites, elles sont conditionnellement stables. Pour ne pas être limité par le pas de temps de la convection, il est possible de faire la convection en plusieurs pas (sous pas de convection). L'étape de convection peut être configurée comme suit:

Méthode #1:

ale.getTransferStep().setSubStepping(n)

divise chaque pas de convection en n pas.

Méthode #2:

ale.getTransferStep().setAutoSubStepping(cflMax)
ale.getTransferStep().setCriticalCfl(cflCrit)

divise automatiquement le pas de convection pour que le cflMax (condition de Courant) ne soit jamais atteint. Pour info le CFL doit toujours être inférieur à 1 pour être stable. On peut donc spécifier ici un facteur de sécurité. Par exemple 0.8 ou 0.5. La valeur de cflCrit est la valeur pour laquelle on suppose que la convection est beaucoup trop importante et on divisera alors le pas de temps si cette valeur est atteinte (ce qui revient à rejeter le pas Lagrangien même s'il a convergé et redémarrer avec un incrément de charge plus faible).

Lorsque la matière est susceptible d'entrer dans le maillage, il est parfois nécessaire de définir des conditions aux limites. Par défaut, les valeurs des conditions aux limites sont égales aux valeurs des grandeurs sur la frontière. En d'autres mots, si de la matière entre dans le maillage par une frontière, elle aura les même caractéristique que cette frontière.

Si on veut appliquer d'autres valeurs (par exemple forcer de la matière vierge en entrée), il faut alors explicitement définir de nouveaux éléments de conditions aux limites via une interaction dédiée.

Cette interaction est de type BcInteraction et les éléments de condition aux limites s'appellent BoundaryCell:

prp = ElementProperties(BoundaryCell)   # propriétés
prp.put(MATERIAL, 1)                    # matériau (le même que celui de la région)

bci = BcInteraction(1)                  # le générateur de ''BoundaryCell''
bci.push(curveset(5))                   # frontière concernée (ici courbe #5)
bci.addProperty(prp)                    # lien propriété-interaction

region.getInteractionSet().add(bci)     # ajout au niveau de la région ALE concernée.

Ces lignes vont définir des éléments de condition aux limites sur la ligne #5. Ceux-ci se comporteront par défaut exactement comme s'ils n'existaient pas.

Pour profiter de ces éléments, il faut ensuite spécifier des valeurs. Ceci se fait à l'aide d'un LoadingSet au niveau de la région ALE. Par exemple:

region.getLoadingSet().define(bci, IF_FTOTAL_XY, 0.0)
region.getLoadingSet().define(bci, IF_FTOTAL_YX, 0.0)

Ces commandes imposent une valeur 0.0 en entrée aux composantes de cisaillement du tenseur gradient de défos.

Après convection, certaines opérations doivent être effectuées. Il s'agit par exemple de recalculer le glissement des contacts avec frottement ou recalculer une matrice des masses sur le nouveau maillage.

Chaque region execute une commande de post-processing. Pour ce qui est des contacts, il faut explicitement spécifier que l'on veut les rééquilibrer. Cela se fait par la commande:

ale.getPostStep().add(ci1)

où ci1 est une interaction de contact.